C++ FEInterpolation类代码示例

void MixedGradientPressureWeakPeriodic :: computeStress(FloatArray &sigmaDev, FloatArray &tractions, double rve_size)
{
    FloatMatrix mMatrix;
    FloatArray normal, coords, t;

    int nsd = domain->giveNumberOfSpatialDimensions();
    Set *set = this->giveDomain()->giveSet(this->set);
    const IntArray &boundaries = set->giveBoundaryList();
    // Reminder: sigma = int t * n dA, where t = sum( C_i * n t_i ).
    // This loop will construct sigma in matrix form.

    FloatMatrix sigma;

    for ( int pos = 1; pos <= boundaries.giveSize() / 2; ++pos ) {
        Element *el = this->giveDomain()->giveElement( boundaries.at(pos * 2 - 1) );
        int boundary = boundaries.at(pos * 2);

        FEInterpolation *interp = el->giveInterpolation(); // Geometry interpolation. The displacements or velocities must have the same interpolation scheme (on the boundary at least).

        int maxorder = this->order + interp->giveInterpolationOrder() * 3;
        std :: unique_ptr< IntegrationRule >ir( interp->giveBoundaryIntegrationRule(maxorder, boundary) );

        for ( GaussPoint *gp: *ir ) {
            const FloatArray &lcoords = gp->giveNaturalCoordinates();
            FEIElementGeometryWrapper cellgeo(el);

            double detJ = interp->boundaryEvalNormal(normal, boundary, lcoords, cellgeo);
            // Compute v_m = d_dev . x
            interp->boundaryLocal2Global(coords, boundary, lcoords, cellgeo);

            this->constructMMatrix(mMatrix, coords, normal);
            t.beProductOf(mMatrix, tractions);
            sigma.plusDyadUnsym(t, coords, detJ * gp->giveWeight());
        }
    }
    sigma.times(1. / rve_size);

    double pressure = 0.;
    for ( int i = 1; i <= nsd; i++ ) {
        pressure += sigma.at(i, i);
    }
    pressure /= 3; // Not 100% sure about this for 2D cases.
    if ( nsd == 3 ) {
        sigmaDev.resize(6);
        sigmaDev.at(1) = sigma.at(1, 1) - pressure;
        sigmaDev.at(2) = sigma.at(2, 2) - pressure;
        sigmaDev.at(3) = sigma.at(3, 3) - pressure;
        sigmaDev.at(4) = 0.5 * ( sigma.at(2, 3) + sigma.at(3, 2) );
        sigmaDev.at(5) = 0.5 * ( sigma.at(1, 3) + sigma.at(3, 1) );
        sigmaDev.at(6) = 0.5 * ( sigma.at(1, 2) + sigma.at(2, 1) );
    } else if ( nsd == 2 ) {
        sigmaDev.resize(3);
        sigmaDev.at(1) = sigma.at(1, 1) - pressure;
        sigmaDev.at(2) = sigma.at(2, 2) - pressure;
        sigmaDev.at(3) = 0.5 * ( sigma.at(1, 2) + sigma.at(2, 1) );
    } else {
        sigmaDev.resize(1);
        sigmaDev.at(1) = sigma.at(1, 1) - pressure;
    }
}

void MixedGradientPressureWeakPeriodic :: integrateTractionXTangent(FloatMatrix &answer, Element *el, int boundary)
{
    // Computes the integral: int dt . dx_m dA
    FloatMatrix mMatrix;
    FloatArray normal, coords, vM_vol;

    FEInterpolation *interp = el->giveInterpolation(); // Geometry interpolation. The displacements or velocities must have the same interpolation scheme (on the boundary at least).

    int maxorder = this->order + interp->giveInterpolationOrder() * 3;
    std :: unique_ptr< IntegrationRule >ir( interp->giveBoundaryIntegrationRule(maxorder, boundary) );

    FloatArray tmpAnswer;
    for ( GaussPoint *gp: *ir ) {
        const FloatArray &lcoords = gp->giveNaturalCoordinates();
        FEIElementGeometryWrapper cellgeo(el);

        double detJ = interp->boundaryEvalNormal(normal, boundary, lcoords, cellgeo);
        interp->boundaryLocal2Global(coords, boundary, lcoords, cellgeo);

        vM_vol.beScaled(1.0/3.0, coords);
        this->constructMMatrix(mMatrix, coords, normal);

        tmpAnswer.plusProduct(mMatrix, vM_vol, detJ * gp->giveWeight());
    }
    answer.resize(tmpAnswer.giveSize(), 1);
    answer.setColumn(tmpAnswer, 1);
}

void MixedGradientPressureWeakPeriodic :: integrateTractionDev(FloatArray &answer, Element *el, int boundary, const FloatMatrix &ddev)
{
    // Computes the integral: int dt . dx dA
    FloatMatrix mMatrix;
    FloatArray normal, coords, vM_dev;

    FEInterpolation *interp = el->giveInterpolation(); // Geometry interpolation. The displacements or velocities must have the same interpolation scheme (on the boundary at least).

    int maxorder = this->order + interp->giveInterpolationOrder() * 3;
    std :: unique_ptr< IntegrationRule >ir( interp->giveBoundaryIntegrationRule(maxorder, boundary) );
    answer.clear();

    for ( GaussPoint *gp: *ir ) {
        const FloatArray &lcoords = gp->giveNaturalCoordinates();
        FEIElementGeometryWrapper cellgeo(el);

        double detJ = interp->boundaryEvalNormal(normal, boundary, lcoords, cellgeo);
        // Compute v_m = d_dev . x
        interp->boundaryLocal2Global(coords, boundary, lcoords, cellgeo);
        vM_dev.beProductOf(ddev, coords);

        this->constructMMatrix(mMatrix, coords, normal);

        answer.plusProduct(mMatrix, vM_dev, detJ * gp->giveWeight());
    }
}

int
DofManValueField :: evaluateAt(FloatArray &answer, const FloatArray &coords, ValueModeType mode, TimeStep *tStep)
{
    int result = 0; // assume ok
    FloatArray lc, n;

    // request element containing target point
    Element *elem = this->domain->giveSpatialLocalizer()->giveElementContainingPoint(coords);
    if ( elem ) { // ok element containing target point found
        FEInterpolation *interp = elem->giveInterpolation();
        if ( interp ) {
            // map target point to element local coordinates
            if ( interp->global2local( lc, coords, FEIElementGeometryWrapper(elem) ) ) {
                // evaluate interpolation functions at target point
                interp->evalN( n, lc, FEIElementGeometryWrapper(elem) );
                // loop over element nodes
                for ( int i = 1; i <= n.giveSize(); i++ ) {
                    // multiply nodal value by value of corresponding shape function and add this to answer
                    answer.add( n.at(i), this->dmanvallist[elem->giveDofManagerNumber(i)-1] );
                }
            } else { // mapping from global to local coordinates failed
                result = 1; // failed
            }
        } else {  // element without interpolation
            result = 1; // failed
        }
    } else { // no element containing given point found
        result = 1; // failed
    }
    return result;
}

void
AxisymElement :: computeBmatrixAt(GaussPoint *gp, FloatMatrix &answer, int li, int ui)
//
// Returns the [ 6 x (nno*2) ] strain-displacement matrix {B} of the receiver,
// evaluated at gp.
// (epsilon_x,epsilon_y,...,Gamma_xy) = B . r
// r = ( u1,v1,u2,v2,u3,v3,u4,v4)
{

    FEInterpolation *interp = this->giveInterpolation();
     
    FloatArray N;
    interp->evalN( N, gp->giveNaturalCoordinates(), *this->giveCellGeometryWrapper() );
    double r = 0.0;
    for ( int i = 1; i <= this->giveNumberOfDofManagers(); i++ ) {
        double x = this->giveNode(i)->giveCoordinate(1);
        r += x * N.at(i);
    } 
    
    FloatMatrix dNdx;
    interp->evaldNdx( dNdx, gp->giveNaturalCoordinates(), *this->giveCellGeometryWrapper() );
    answer.resize(6, dNdx.giveNumberOfRows() * 2);
    answer.zero();

    for ( int i = 1; i <= dNdx.giveNumberOfRows(); i++ ) {
        answer.at(1, i * 2 - 1) = dNdx.at(i, 1);
        answer.at(2, i * 2 - 0) = dNdx.at(i, 2);
        answer.at(3, i * 2 - 1) = N.at(i) / r;
        answer.at(6, 2 * i - 1) = dNdx.at(i, 2);
        answer.at(6, 2 * i - 0) = dNdx.at(i, 1);
    }
  
}

void Space3dStructuralElementEvaluator :: computeBMatrixAt(FloatMatrix &answer, GaussPoint *gp)
{
    FloatMatrix d;
    Element *element = this->giveElement();
    FEInterpolation *interp = element->giveInterpolation();
    // this uses FEInterpolation::nodes2coords - quite inefficient in this case (large num of dofmans)
    interp->evaldNdx( d, gp->giveNaturalCoordinates(), FEIIGAElementGeometryWrapper( element, gp->giveIntegrationRule()->giveKnotSpan() ) );


    answer.resize(6, d.giveNumberOfRows() * 3);
    answer.zero();

    for ( int i = 1; i <= d.giveNumberOfRows(); i++ ) {
        answer.at(1, i * 3 - 2) = d.at(i, 1);
        answer.at(2, i * 3 - 1) = d.at(i, 2);
        answer.at(3, i * 3 - 0) = d.at(i, 3);

        answer.at(4, 3 * i - 1) = d.at(i, 3);
        answer.at(4, 3 * i - 0) = d.at(i, 2);

        answer.at(5, 3 * i - 2) = d.at(i, 3);
        answer.at(5, 3 * i - 0) = d.at(i, 1);

        answer.at(6, 3 * i - 2) = d.at(i, 2);
        answer.at(6, 3 * i - 1) = d.at(i, 1);
    }
}

bool EnrichmentItem :: tipIsTouchingEI(const TipInfo &iTipInfo)
{
    double tol = 1.0e-9;
    SpatialLocalizer *localizer = giveDomain()->giveSpatialLocalizer();

    Element *tipEl = localizer->giveElementContainingPoint(iTipInfo.mGlobalCoord);
    if ( tipEl != NULL ) {
        // Check if the candidate tip is located on the current crack
        FloatArray N;
        FloatArray locCoord;
        tipEl->computeLocalCoordinates(locCoord, iTipInfo.mGlobalCoord);
        FEInterpolation *interp = tipEl->giveInterpolation();
        interp->evalN( N, locCoord, FEIElementGeometryWrapper(tipEl) );

        double normalSignDist;
        evalLevelSetNormal( normalSignDist, iTipInfo.mGlobalCoord, N, tipEl->giveDofManArray() );

        double tangSignDist;
        evalLevelSetTangential( tangSignDist, iTipInfo.mGlobalCoord, N, tipEl->giveDofManArray() );

        if ( fabs(normalSignDist) < tol && tangSignDist > tol ) {
            return true;
        }
    }

    return false;
}

void PlaneStressStructuralElementEvaluator :: computeNMatrixAt(FloatMatrix &answer, GaussPoint *gp)
{
    FloatArray N;
    FEInterpolation *interp = gp->giveElement()->giveInterpolation();
    interp->evalN( N, gp->giveNaturalCoordinates(), FEIIGAElementGeometryWrapper( gp->giveElement(), gp->giveIntegrationRule()->giveKnotSpan() ) );
    answer.beNMatrixOf(N, 2);
}

void
Structural3DElement :: computeBHmatrixAt(GaussPoint *gp, FloatMatrix &answer)
// Returns the [ 9 x (nno * 3) ] displacement gradient matrix {BH} of the receiver,
// evaluated at gp.
// BH matrix  -  9 rows : du/dx, dv/dy, dw/dz, dv/dz, du/dz, du/dy, dw/dy, dw/dx, dv/dx
{
    FEInterpolation *interp = this->giveInterpolation();
    FloatMatrix dNdx; 
    interp->evaldNdx( dNdx, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );
    
    answer.resize(9, dNdx.giveNumberOfRows() * 3);
    answer.zero();

    for ( int i = 1; i <= dNdx.giveNumberOfRows(); i++ ) {
        answer.at(1, 3 * i - 2) = dNdx.at(i, 1);     // du/dx
        answer.at(2, 3 * i - 1) = dNdx.at(i, 2);     // dv/dy
        answer.at(3, 3 * i - 0) = dNdx.at(i, 3);     // dw/dz
        answer.at(4, 3 * i - 1) = dNdx.at(i, 3);     // dv/dz
        answer.at(7, 3 * i - 0) = dNdx.at(i, 2);     // dw/dy
        answer.at(5, 3 * i - 2) = dNdx.at(i, 3);     // du/dz
        answer.at(8, 3 * i - 0) = dNdx.at(i, 1);     // dw/dx
        answer.at(6, 3 * i - 2) = dNdx.at(i, 2);     // du/dy
        answer.at(9, 3 * i - 1) = dNdx.at(i, 1);     // dv/dx
    }

}

void
MatlabExportModule :: doOutputHomogenizeDofIDs(TimeStep *tStep,    FILE *FID)
{

    std :: vector <FloatArray*> HomQuantities;
    double Vol = 0.0;

    // Initialize vector of arrays constaining homogenized quantities
    HomQuantities.resize(internalVarsToExport.giveSize());

    for (int j=0; j<internalVarsToExport.giveSize(); j++) {
        HomQuantities.at(j) = new FloatArray;
    }

    int nelem = this->elList.giveSize();
    for (int i = 1; i<=nelem; i++) {
        Element *e = this->emodel->giveDomain(1)->giveElement(elList.at(i));
        FEInterpolation *Interpolation = e->giveInterpolation();

        Vol = Vol + e->computeVolumeAreaOrLength();

        for ( GaussPoint *gp: *e->giveDefaultIntegrationRulePtr() ) {

            for (int j=0; j<internalVarsToExport.giveSize(); j++) {
                FloatArray elementValues;
                e->giveIPValue(elementValues, gp, (InternalStateType) internalVarsToExport(j), tStep);
                double detJ=fabs(Interpolation->giveTransformationJacobian( gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(e)));

                elementValues.times(gp->giveWeight()*detJ);
                if (HomQuantities.at(j)->giveSize() == 0) {
                    HomQuantities.at(j)->resize(elementValues.giveSize());
                    HomQuantities.at(j)->zero();
                };
                HomQuantities.at(j)->add(elementValues);
            }
        }
    }


    if (noscaling) Vol=1.0;

    for ( std :: size_t i = 0; i < HomQuantities.size(); i ++) {
        FloatArray *thisIS;
        thisIS = HomQuantities.at(i);
        thisIS->times(1.0/Vol);
        fprintf(FID, "\tspecials.%s = [", __InternalStateTypeToString ( InternalStateType (internalVarsToExport(i)) ) );

        for (int j = 0; j<thisIS->giveSize(); j++) {
            fprintf(FID, "%e", thisIS->at(j+1));
            if (j!=(thisIS->giveSize()-1) ) {
                fprintf(FID, ", ");
            }
        }
        fprintf(FID, "];\n");
        delete HomQuantities.at(i);
    }

}

/* 3D Space Elements */
void Space3dStructuralElementEvaluator :: computeNMatrixAt(FloatMatrix &answer, GaussPoint *gp)
{
    FloatArray N;
    Element *element = this->giveElement();
    FEInterpolation *interp = element->giveInterpolation();

    interp->evalN( N, gp->giveNaturalCoordinates(), FEIIGAElementGeometryWrapper( element, gp->giveIntegrationRule()->giveKnotSpan() ) );

    answer.beNMatrixOf(N, 3);
}

void
PrescribedMean :: computeDomainSize()
{
    if (domainSize > 0.0) return;


    if (elementEdges) {
        IntArray setList = ((GeneralBoundaryCondition *)this)->giveDomain()->giveSet(set)->giveBoundaryList();

        elements.resize(setList.giveSize() / 2);
        sides.resize(setList.giveSize() / 2);

        for (int i=1; i<=setList.giveSize(); i=i+2) {
            elements.at(i/2+1) = setList.at(i);
            sides.at(i/2+1) = setList.at(i+1);
        }
    } else {
        IntArray setList = ((GeneralBoundaryCondition *)this)->giveDomain()->giveSet(set)->giveElementList();
        elements = setList;
    }

    domainSize = 0.0;

    for ( int i=1; i<=elements.giveSize(); i++ ) {
        int elementID = elements.at(i);
        Element *thisElement = this->giveDomain()->giveElement(elementID);
        FEInterpolation *interpolator = thisElement->giveInterpolation(DofIDItem(dofid));

        IntegrationRule *iRule;

        if (elementEdges) {
            iRule = interpolator->giveBoundaryIntegrationRule(3, sides.at(i));
        } else {
            iRule = interpolator->giveIntegrationRule(3);
        }

        for ( GaussPoint * gp: * iRule ) {
            FloatArray lcoords = gp->giveNaturalCoordinates();

            double detJ;
            if (elementEdges) {
                detJ = fabs ( interpolator->boundaryGiveTransformationJacobian(sides.at(i), lcoords, FEIElementGeometryWrapper(thisElement)) );
            } else {
                detJ = fabs ( interpolator->giveTransformationJacobian(lcoords, FEIElementGeometryWrapper(thisElement)) );
            }
            domainSize = domainSize + detJ*gp->giveWeight();
        }

        delete iRule;
    }

    printf("%f\n", domainSize);

}

double MixedGradientPressureBC :: domainSize()
{
    int nsd = this->domain->giveNumberOfSpatialDimensions();
    double domain_size = 0.0;
    // This requires the boundary to be consistent and ordered correctly.
    Set *set = this->giveDomain()->giveSet(this->set);
    const IntArray &boundaries = set->giveBoundaryList();

    for ( int pos = 1; pos <= boundaries.giveSize() / 2; ++pos ) {
        Element *e = this->giveDomain()->giveElement( boundaries.at(pos * 2 - 1) );
        int boundary = boundaries.at(pos * 2);
        FEInterpolation *fei = e->giveInterpolation();
        domain_size += fei->evalNXIntegral( boundary, FEIElementGeometryWrapper(e) );
    }
    return domain_size / nsd;
}

void
IntElLine1PhF :: computeCovarBaseVectorAt(IntegrationPoint *ip, FloatArray &G)
{
    FloatMatrix dNdxi;
    FEInterpolation *interp = this->giveInterpolation();
    interp->evaldNdxi( dNdxi, ip->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );
    G.resize(2);
    G.zero();
    int numNodes = this->giveNumberOfNodes();
    for ( int i = 1; i <= dNdxi.giveNumberOfRows(); i++ ) {
        double X1_i = 0.5 * ( this->giveNode(i)->giveCoordinate(1) + this->giveNode(i + numNodes / 2)->giveCoordinate(1) ); // (mean) point on the fictious mid surface
        double X2_i = 0.5 * ( this->giveNode(i)->giveCoordinate(2) + this->giveNode(i + numNodes / 2)->giveCoordinate(2) );
        G.at(1) += dNdxi.at(i, 1) * X1_i;
        G.at(2) += dNdxi.at(i, 1) * X2_i;
    }
}

double TransportGradientPeriodic :: domainSize(Domain *d, int setNum)
{
    int nsd = d->giveNumberOfSpatialDimensions();
    double domain_size = 0.0;
    // This requires the boundary to be consistent and ordered correctly.
    Set *set = d->giveSet(setNum);
    const IntArray &boundaries = set->giveBoundaryList();

    for ( int pos = 1; pos <= boundaries.giveSize() / 2; ++pos ) {
        Element *e = d->giveElement( boundaries.at(pos * 2 - 1) );
        int boundary = boundaries.at(pos * 2);
        FEInterpolation *fei = e->giveInterpolation();
        domain_size += fei->evalNXIntegral( boundary, FEIElementGeometryWrapper(e) );
    }
    return fabs(domain_size / nsd);
}