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

double
FEI3dLineLin :: evaldNdx(FloatMatrix &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
{
    ///@todo Not clear what this function should return. Just dNds would make sense if the caller defines a local coordinate system.
    FloatArray vec;
    vec.beDifferenceOf(*cellgeo.giveVertexCoordinates(2), *cellgeo.giveVertexCoordinates(1));

    double detJ = vec.computeSquaredNorm() * 0.5;
    double l2_inv = 0.5 / detJ;
    answer.resize(2, 3);

    answer.at(1, 1) = -vec.at(1)*l2_inv;
    answer.at(2, 1) =  vec.at(1)*l2_inv;
    answer.at(1, 2) = -vec.at(2)*l2_inv;
    answer.at(2, 2) =  vec.at(2)*l2_inv;
    answer.at(1, 3) = -vec.at(3)*l2_inv;
    answer.at(2, 3) =  vec.at(3)*l2_inv;

    return detJ;
}

bool XfemStructuralElementInterface :: XfemElementInterface_updateIntegrationRule()
{
    const double tol2 = 1.0e-18;

    bool partitionSucceeded = false;


    if ( mpCZMat != NULL ) {
        mpCZIntegrationRules.clear();
        mCZEnrItemIndices.clear();
        mCZTouchingEnrItemIndices.clear();
    }

    XfemManager *xMan = this->element->giveDomain()->giveXfemManager();
    if ( xMan->isElementEnriched(element) ) {
        if ( mpCZMat == NULL && mCZMaterialNum > 0 ) {
            initializeCZMaterial();
        }


        MaterialMode matMode = element->giveMaterialMode();

        bool firstIntersection = true;

        std :: vector< std :: vector< FloatArray > >pointPartitions;
        mSubTri.clear();

        std :: vector< int >enrichingEIs;
        int elPlaceInArray = xMan->giveDomain()->giveElementPlaceInArray( element->giveGlobalNumber() );
        xMan->giveElementEnrichmentItemIndices(enrichingEIs, elPlaceInArray);


        for ( size_t p = 0; p < enrichingEIs.size(); p++ ) {
            // Index of current ei
            int eiIndex = enrichingEIs [ p ];

            // Indices of other ei interaction with this ei through intersection enrichment fronts.
            std :: vector< int >touchingEiIndices;
            giveIntersectionsTouchingCrack(touchingEiIndices, enrichingEIs, eiIndex, * xMan);

            if ( firstIntersection ) {
                // Get the points describing each subdivision of the element
                double startXi, endXi;
                bool intersection = false;
                this->XfemElementInterface_prepareNodesForDelaunay(pointPartitions, startXi, endXi, eiIndex, intersection);

                if ( intersection ) {
                    firstIntersection = false;

                    // Use XfemElementInterface_partitionElement to subdivide the element
                    for ( int i = 0; i < int ( pointPartitions.size() ); i++ ) {
                        // Triangulate the subdivisions
                        this->XfemElementInterface_partitionElement(mSubTri, pointPartitions [ i ]);
                    }


                    if ( mpCZMat != NULL ) {
                        Crack *crack = dynamic_cast< Crack * >( xMan->giveEnrichmentItem(eiIndex) );
                        if ( crack == NULL ) {
                            OOFEM_ERROR("Cohesive zones are only available for cracks.")
                        }

                        // We have xi_s and xi_e. Fetch sub polygon.
                        std :: vector< FloatArray >crackPolygon;
                        crack->giveSubPolygon(crackPolygon, startXi, endXi);

                        ///////////////////////////////////
                        // Add cohesive zone Gauss points
                        size_t numSeg = crackPolygon.size() - 1;

                        for ( size_t segIndex = 0; segIndex < numSeg; segIndex++ ) {
                            int czRuleNum = 1;
                            mpCZIntegrationRules.emplace_back( new GaussIntegrationRule(czRuleNum, element) );

                            // Add index of current ei
                            mCZEnrItemIndices.push_back(eiIndex);

                            // Add indices of other ei, that cause interaction through
                            // intersection enrichment fronts
                            mCZTouchingEnrItemIndices.push_back(touchingEiIndices);

                            // Compute crack normal
                            FloatArray crackTang;
                            crackTang.beDifferenceOf(crackPolygon [ segIndex + 1 ], crackPolygon [ segIndex ]);

                            if ( crackTang.computeSquaredNorm() > tol2 ) {
                                crackTang.normalize();
                            }

                            FloatArray crackNormal = {
                                -crackTang.at(2), crackTang.at(1)
                            };

                            mpCZIntegrationRules [ segIndex ]->SetUpPointsOn2DEmbeddedLine(mCSNumGaussPoints, matMode,
                                                                                           crackPolygon [ segIndex ], crackPolygon [ segIndex + 1 ]);

                            for ( GaussPoint *gp: *mpCZIntegrationRules [ segIndex ] ) {
                                double gw = gp->giveWeight();
                                double segLength = crackPolygon [ segIndex ].distance(crackPolygon [ segIndex + 1 ]);
                                gw *= 0.5 * segLength;
//.........这里部分代码省略.........

void PolygonLine :: computeNormalSignDist(double &oDist, const FloatArray &iPoint) const
{
    const FloatArray &point = {iPoint[0], iPoint[1]};

    oDist = std :: numeric_limits< double > :: max();
    int numSeg = this->giveNrVertices() - 1;

    // TODO: This can probably be done in a nicer way.
    // Ensure that we work in 2d.
    const int dim = 2;

    for ( int segId = 1; segId <= numSeg; segId++ ) {
        // Crack segment
        const FloatArray &crackP1( this->giveVertex ( segId ) );

        const FloatArray &crackP2( this->giveVertex ( segId + 1 ) );

        double dist2 = 0.0;
        if ( segId == 1 ) {
            // Vector from start P1 to point X
            FloatArray u = {point.at(1) - crackP1.at(1), point.at(2) - crackP1.at(2)};

            // Line tangent vector
            FloatArray t = {crackP2.at(1) - crackP1.at(1), crackP2.at(2) - crackP1.at(2)};
            double l2 = t.computeSquaredNorm();

            if ( l2 > 0.0 ) {
                double l = t.normalize();
                double s = dot(u, t);

                if ( s > l ) {
                    // X is closest to P2
                    dist2 = point.distance_square(crackP2);
                } else {
                    double xi = s / l;
                    FloatArray q = ( 1.0 - xi ) * crackP1 + xi * crackP2;
                    dist2 = point.distance_square(q);
                }
            } else {
                // If the points P1 and P2 coincide,
                // we can compute the distance to any
                // of these points.
                dist2 = point.distance_square(crackP1);
            }
        } else if ( segId == numSeg ) {
            // Vector from start P1 to point X
            FloatArray u = {point.at(1) - crackP1.at(1), point.at(2) - crackP1.at(2)};

            // Line tangent vector
            FloatArray t = {crackP2.at(1) - crackP1.at(1), crackP2.at(2) - crackP1.at(2)};
            double l2 = t.computeSquaredNorm();

            if ( l2 > 0.0 ) {
                double l = t.normalize();
                double s = dot(u, t);

                if ( s < 0.0 ) {
                    // X is closest to P1
                    dist2 = point.distance_square(crackP1);
                } else {
                    double xi = s / l;
                    FloatArray q = ( 1.0 - xi ) * crackP1 + xi * crackP2;
                    dist2 = point.distance_square(q);
                }
            } else {
                // If the points P1 and P2 coincide,
                // we can compute the distance to any
                // of these points.
                dist2 = point.distance_square(crackP1);
            }
        } else {
            double arcPos = -1.0, dummy;
            dist2 = point.distance_square(crackP1, crackP2, arcPos, dummy);
        }

        if ( dist2 < oDist*oDist ) {
            FloatArray lineToP;
            lineToP.beDifferenceOf(point, crackP1, dim);

            FloatArray t;
            t.beDifferenceOf(crackP2, crackP1, dim);

            FloatArray n = {-t.at(2), t.at(1)};

            oDist = sgn( lineToP.dotProduct(n) ) * sqrt(dist2);
        }
    }
}

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

void
StructuralEngngModel :: giveInternalForces(FloatArray &answer, bool normFlag, int di, TimeStep *stepN)
{
    // Simply assembles contributions from each element in domain
    Element *element;
    IntArray loc;
    FloatArray charVec;
    FloatMatrix R;
    int nelems;
    EModelDefaultEquationNumbering dn;

    answer.resize( this->giveNumberOfEquations( EID_MomentumBalance ) );
    answer.zero();
    if ( normFlag ) internalForcesEBENorm = 0.0;

    // Update solution state counter
    stepN->incrementStateCounter();

#ifdef __PARALLEL_MODE
    if ( isParallel() ) {
        exchangeRemoteElementData( RemoteElementExchangeTag  );
    }
#endif

    nelems = this->giveDomain(di)->giveNumberOfElements();
    this->timer.resumeTimer(EngngModelTimer :: EMTT_NetComputationalStepTimer);

    for ( int i = 1; i <= nelems; i++ ) {
        element = this->giveDomain(di)->giveElement(i);
#ifdef __PARALLEL_MODE
        // Skip remote elements (these are used as mirrors of remote elements on other domains
        // when nonlocal constitutive models are used. Their introduction is necessary to
        // allow local averaging on domains without fine grain communication between domains).
        if ( element->giveParallelMode() == Element_remote ) continue;

#endif
        element->giveLocationArray( loc, EID_MomentumBalance, dn );
        element->giveCharacteristicVector( charVec, InternalForcesVector, VM_Total, stepN );
        if ( element->giveRotationMatrix(R, EID_MomentumBalance) ) {
            charVec.rotatedWith(R, 't');
        }
        answer.assemble(charVec, loc);

        // Compute element norm contribution.
        if ( normFlag ) internalForcesEBENorm += charVec.computeSquaredNorm();
    }

    this->timer.pauseTimer( EngngModelTimer :: EMTT_NetComputationalStepTimer );

#ifdef __PARALLEL_MODE
    this->updateSharedDofManagers(answer, InternalForcesExchangeTag);

    if ( normFlag ) {
        // Exchange norm contributions.
        double localEBENorm = internalForcesEBENorm;
        internalForcesEBENorm = 0.0;
        MPI_Allreduce(& localEBENorm, & internalForcesEBENorm, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
    }
#endif

    // Remember last internal vars update time stamp.
    internalVarUpdateStamp = stepN->giveSolutionStateCounter();
}

//.........这里部分代码省略.........
            
            OOFEM_LOG_INFO( "  %s:", __DofIDItemToString((DofIDItem)dg).c_str() );

            if ( rtolf.at(1) > 0.0 ) {
                //  compute a relative error norm
                if ( ( dg_totalLoadLevel.at(dg) + internalForcesEBENorm.at(dg) ) > nrsolver_ERROR_NORM_SMALL_NUM ) {
                    forceErr = sqrt( dg_forceErr.at(dg) / ( dg_totalLoadLevel.at(dg) + internalForcesEBENorm.at(dg) ) );
                } else {
                    // If both external forces and internal ebe norms are zero, then the residual must be zero.
                    //zeroNorm = true; // Warning about this afterwards.
                    zeroFNorm = true;
                    forceErr = sqrt( dg_forceErr.at(dg) );
                }

                if ( forceErr > rtolf.at(1) * NRSOLVER_MAX_REL_ERROR_BOUND ) {
                    errorOutOfRange = true;
                }
                if ( forceErr > rtolf.at(1) ) {
                    answer = false;
                }
                OOFEM_LOG_INFO( zeroFNorm ? " *%.3e" : "  %.3e", forceErr );
            }

            if ( rtold.at(1) > 0.0 ) {
                // compute displacement error
                if ( dg_totalDisp.at(dg) >  nrsolver_ERROR_NORM_SMALL_NUM ) {
                    dispErr = sqrt( dg_dispErr.at(dg) / dg_totalDisp.at(dg) );
                } else {
                    ///@todo This is almost always the case for displacement error. nrsolveR_ERROR_NORM_SMALL_NUM is no good.
                    //zeroNorm = true; // Warning about this afterwards.
                    //zeroDNorm = true;
                    dispErr = sqrt( dg_dispErr.at(dg) );
                }
                if ( dispErr  > rtold.at(1) * NRSOLVER_MAX_REL_ERROR_BOUND ) {
                    errorOutOfRange = true;
                }
                if ( dispErr > rtold.at(1) ) {
                    answer = false;
                }
                OOFEM_LOG_INFO( zeroDNorm ? " *%.3e" : "  %.3e", dispErr );
            }
        }
        OOFEM_LOG_INFO("\n");
        //if ( zeroNorm ) OOFEM_WARNING("NRSolver :: checkConvergence - Had to resort to absolute error measure (marked by *)");
    } else { // No dof grouping
        double dXX, dXdX;
        
        if ( engngModel->giveProblemScale() == macroScale ) {
            OOFEM_LOG_INFO("NRSolver:     %-15d", nite);
        } else {
            OOFEM_LOG_INFO("  NRSolver:     %-15d", nite);
        }

 #ifdef __PARALLEL_MODE
        forceErr = parallel_context->norm(rhs); forceErr *= forceErr;
        dXX = parallel_context->localNorm(X); dXX *= dXX; // Note: Solutions are always total global values (natural distribution makes little sense for the solution)
        dXdX = parallel_context->localNorm(ddX); dXdX *= dXdX;
 #else
        forceErr = rhs.computeSquaredNorm();
        dXX = X.computeSquaredNorm();
        dXdX = ddX.computeSquaredNorm();
 #endif
        if ( rtolf.at(1) > 0.0 ) {
            // we compute a relative error norm
            if ( ( RRT + internalForcesEBENorm.at(1) ) > nrsolver_ERROR_NORM_SMALL_NUM ) {
                forceErr = sqrt( forceErr / ( RRT + internalForcesEBENorm.at(1) ) );
            } else {
                forceErr = sqrt( forceErr ); // absolute norm as last resort
            }
            if ( fabs(forceErr) > rtolf.at(1) * NRSOLVER_MAX_REL_ERROR_BOUND ) {
                errorOutOfRange = true;
            }
            if ( fabs(forceErr) > rtolf.at(1) ) {
                answer = false;
            }
            OOFEM_LOG_INFO(" %-15e", forceErr);
        }

        if ( rtold.at(1) > 0.0 ) {
            // compute displacement error
            // err is relative displacement change
            if ( dXX > nrsolver_ERROR_NORM_SMALL_NUM ) {
                dispErr = sqrt( dXdX / dXX );
            } else {
                dispErr = sqrt( dXdX );
            }
            if ( fabs(dispErr)  > rtold.at(1) * NRSOLVER_MAX_REL_ERROR_BOUND ) {
                errorOutOfRange = true;
            }
            if ( fabs(dispErr)  > rtold.at(1) ) {
                answer = false;
            }
            OOFEM_LOG_INFO(" %-15e", dispErr);
        }

        OOFEM_LOG_INFO("\n");
    } // end default case (all dofs contributing)

    return answer;
}