C++ DofManager::giveCoordinate方法代码示例

double
UserDefDirichletBC :: give(Dof *dof, ValueModeType mode, TimeStep *stepN)
{
    double factor = this->giveLoadTimeFunction()->evaluate(stepN, mode);
    DofManager *dMan = dof->giveDofManager();


    /*
     * The Python function takes two input arguments:
     * 	1) An array with node coordinates
     * 	2) The dof id
     */
    int numArgs = 3;

    // Create array with node coordinates
    int dim = dMan->giveCoordinates()->giveSize();
    PyObject *pArgArray = PyList_New(dim);

    PyObject *pArgs = PyTuple_New(numArgs);

    for (int i = 0; i < dim; i++) {
        PyList_SET_ITEM(pArgArray, i, PyFloat_FromDouble( dMan->giveCoordinate(i+1) ));
    }
    // PyTuple_SetItem takes over responsibility for objects passed
    // to it -> no DECREF
    PyTuple_SetItem(pArgs, 0, pArgArray);


    // Dof number
    PyObject *pValDofNum = PyLong_FromLong(dof->giveDofID());
    PyTuple_SetItem(pArgs, 1, pValDofNum);


    // Time
    PyObject *pTargetTime = PyFloat_FromDouble( stepN->giveTargetTime() );
    PyTuple_SetItem(pArgs, 2, pTargetTime);

    // Value returned from the Python function
    PyObject *pRetVal;

    if ( PyCallable_Check(mpFunc) ) {
        pRetVal = PyObject_CallObject(mpFunc, pArgs);
    } else {
        OOFEM_ERROR("UserDefDirichletBC :: give: Python function is not callable.");
    }

    // Get return value
    double retVal = 0.0;
    if ( pRetVal != NULL ) {
        retVal = PyFloat_AsDouble(pRetVal);
    }
    else {
        OOFEM_ERROR("UserDefDirichletBC :: give: Failed to fetch Python return value.");
    }

    // Decrement reference count on pointers
    Py_DECREF(pArgs);
    Py_DECREF(pRetVal);

    return retVal*factor;
}

void PLHoopStressCirc :: propagateInterfaces(Domain &iDomain, EnrichmentDomain &ioEnrDom)
{
    // Fetch crack tip data
    TipInfo tipInfoStart, tipInfoEnd;
    ioEnrDom.giveTipInfos(tipInfoStart, tipInfoEnd);
    std :: vector< TipInfo >tipInfo = {tipInfoStart, tipInfoEnd};

    SpatialLocalizer *localizer = iDomain.giveSpatialLocalizer();

    for ( size_t tipIndex = 0; tipIndex < tipInfo.size(); tipIndex++ ) {
        // Construct circle points on an arc from -90 to 90 degrees
        double angle = -90.0 + mAngleInc;
        std :: vector< double >angles;
        while ( angle <= ( 90.0 - mAngleInc ) ) {
            angles.push_back(angle * M_PI / 180.0);
            angle += mAngleInc;
        }

        const FloatArray &xT    = tipInfo [ tipIndex ].mGlobalCoord;
        const FloatArray &t             = tipInfo [ tipIndex ].mTangDir;
        const FloatArray &n             = tipInfo [ tipIndex ].mNormalDir;

        // It is meaningless to propagate a tip that is not inside any element
        Element *el = localizer->giveElementContainingPoint(tipInfo [ tipIndex ].mGlobalCoord);
        if ( el != NULL ) {
            std :: vector< FloatArray >circPoints;

            for ( size_t i = 0; i < angles.size(); i++ ) {
                FloatArray tangent(2);
                tangent.zero();
                tangent.add(cos(angles [ i ]), t);
                tangent.add(sin(angles [ i ]), n);
                tangent.normalize();

                FloatArray x(xT);
                x.add(mRadius, tangent);
                circPoints.push_back(x);
            }



            std :: vector< double >sigTTArray, sigRTArray;

            // Loop over circle points
            for ( size_t pointIndex = 0; pointIndex < circPoints.size(); pointIndex++ ) {
                FloatArray stressVec;

                if ( mUseRadialBasisFunc ) {
                    // Interpolate stress with radial basis functions

                    // Choose a cut-off length l:
                    // take the distance between two nodes in the element containing the
                    // crack tip multiplied by a constant factor.
                    // ( This choice implies that we hope that the element has reasonable
                    // aspect ratio.)
                    const FloatArray &x1 = * ( el->giveDofManager(1)->giveCoordinates() );
                    const FloatArray &x2 = * ( el->giveDofManager(2)->giveCoordinates() );
                    const double l = 1.0 * x1.distance(x2);

                    // Use the octree to get all elements that have
                    // at least one Gauss point in a certain region around the tip.
                    const double searchRadius = 3.0 * l;
                    std :: set< int >elIndices;
                    localizer->giveAllElementsWithIpWithinBox(elIndices, circPoints [ pointIndex ], searchRadius);


                    // Loop over the elements and Gauss points obtained.
                    // Evaluate the interpolation.
                    FloatArray sumQiWiVi;
                    double sumWiVi = 0.0;
                    for ( int elIndex: elIndices ) {
                        Element *gpEl = iDomain.giveElement(elIndex);
                        IntegrationRule *iRule = gpEl->giveDefaultIntegrationRulePtr();

                        for ( GaussPoint *gp_i: *iRule ) {

                            ////////////////////////////////////////
                            // Compute global gp coordinates
                            FloatArray N;
                            FEInterpolation *interp = gpEl->giveInterpolation();
                            interp->evalN( N, * ( gp_i->giveCoordinates() ), FEIElementGeometryWrapper(gpEl) );


                            // Compute global coordinates of Gauss point
                            FloatArray globalCoord(2);
                            globalCoord.zero();

                            for ( int i = 1; i <= gpEl->giveNumberOfDofManagers(); i++ ) {
                                DofManager *dMan = gpEl->giveDofManager(i);
                                globalCoord.at(1) += N.at(i) * dMan->giveCoordinate(1);
                                globalCoord.at(2) += N.at(i) * dMan->giveCoordinate(2);
                            }


                            ////////////////////////////////////////
                            // Compute weight of kernel function

                            FloatArray tipToGP;
                            tipToGP.beDifferenceOf(globalCoord, xT);
                            bool inFrontOfCrack = true;
//.........这里部分代码省略.........

void LSPrimaryVariableMapper :: mapPrimaryVariables(FloatArray &oU, Domain &iOldDom, Domain &iNewDom, ValueModeType iMode, TimeStep &iTStep)
{
    EngngModel *engngMod = iNewDom.giveEngngModel();
    EModelDefaultEquationNumbering num;


    const int dim = iNewDom.giveNumberOfSpatialDimensions();

    int numElNew = iNewDom.giveNumberOfElements();

    // Count dofs
    int numDofsNew = engngMod->giveNumberOfDomainEquations( 1, num );


    oU.resize(numDofsNew);
    oU.zero();

    FloatArray du(numDofsNew);
    du.zero();

    FloatArray res(numDofsNew);

#ifdef __PETSC_MODULE
    PetscSparseMtrx *K = dynamic_cast<PetscSparseMtrx*>( classFactory.createSparseMtrx(SMT_PetscMtrx) );
    SparseLinearSystemNM *solver = classFactory.createSparseLinSolver(ST_Petsc, & iOldDom, engngMod);
#else
    SparseMtrx *K = classFactory.createSparseMtrx(SMT_Skyline);
    SparseLinearSystemNM *solver = classFactory.createSparseLinSolver(ST_Direct, & iOldDom, engngMod);
#endif


    K->buildInternalStructure( engngMod, 1, num );

    int maxIter = 1;

    for ( int iter = 0; iter < maxIter; iter++ ) {
        K->zero();
        res.zero();


        // Contribution from elements
        for ( int elIndex = 1; elIndex <= numElNew; elIndex++ ) {
            StructuralElement *elNew = dynamic_cast< StructuralElement * >( iNewDom.giveElement(elIndex) );
            if ( elNew == NULL ) {
                OOFEM_ERROR("Failed to cast Element new to StructuralElement.");
            }

            ///////////////////////////////////
            // Compute residual

            // Count element dofs
            int numElNodes = elNew->giveNumberOfDofManagers();
            int numElDofs = 0;
            for ( int i = 1; i <= numElNodes; i++ ) {
                numElDofs += elNew->giveDofManager(i)->giveNumberOfDofs();
            }

            FloatArray elRes(numElDofs);
            elRes.zero();

            IntArray elDofsGlob;
            elNew->giveLocationArray( elDofsGlob, num );


            // Loop over Gauss points
            for ( int intRuleInd = 0; intRuleInd < elNew->giveNumberOfIntegrationRules(); intRuleInd++ ) {
                IntegrationRule *iRule = elNew->giveIntegrationRule(intRuleInd);

                for ( GaussPoint *gp: *iRule ) {

                    // New N-matrix
                    FloatMatrix NNew;
                    elNew->computeNmatrixAt(* ( gp->giveNaturalCoordinates() ), NNew);


                    //////////////
                    // Global coordinates of GP
                    const int nDofMan = elNew->giveNumberOfDofManagers();

                    FloatArray Nc;
                    FEInterpolation *interp = elNew->giveInterpolation();
                    const FloatArray &localCoord = * ( gp->giveNaturalCoordinates() );
                    interp->evalN( Nc, localCoord, FEIElementGeometryWrapper(elNew) );

                    const IntArray &elNodes = elNew->giveDofManArray();

                    FloatArray globalCoord(dim);
                    globalCoord.zero();

                    for ( int i = 1; i <= nDofMan; i++ ) {
                        DofManager *dMan = elNew->giveDofManager(i);

                        for ( int j = 1; j <= dim; j++ ) {
                            globalCoord.at(j) += Nc.at(i) * dMan->giveCoordinate(j);
                        }
                    }
                    //////////////


                    // Localize element and point in the old domain
//.........这里部分代码省略.........