C++ SpatialLocalizer::giveElementClosestToPoint方法代码示例

void
FreeWarping :: computeResultAtCenterOfGravity(TimeStep *tStep)
{
    int noCS = this->giveDomain(1)->giveNumberOfCrossSectionModels(); //number of warping Crosssections
    SolutionAtCG.resize(noCS);
    Element *closestElement;
    FloatArray lcoords,  closest, lcg;
    SpatialLocalizer *sp = this->giveDomain(1)->giveSpatialLocalizer();
    sp->init();
    lcoords.resize(2);
    closest.resize(2);
    lcg.resize(2);

    for ( int j = 1; j <= noCS; ++j ) {
        lcg.at(1) = CG.at(j, 1);
        lcg.at(2) = CG.at(j, 2);
        closestElement = sp->giveElementClosestToPoint(lcoords, closest, lcg, 0);

        StructuralElement *sE = dynamic_cast< StructuralElement * >(closestElement);
        FloatArray u, r, b;
        FloatMatrix N;
        sE->computeNmatrixAt(lcoords, N);
        sE->computeVectorOf(VM_Total, tStep, u);
        u.resizeWithValues(3);
        r.beProductOf(N, u);

        SolutionAtCG.at(j) = r.at(1);
    }
}

void
POIExportModule :: exportPrimVarAs(UnknownType valID, FILE *stream, TimeStep *tStep)
{
    Domain *d = emodel->giveDomain(1);
    FloatArray pv, coords(3), lcoords, closest;
    InternalStateValueType type = ISVT_UNDEFINED;

    if ( valID == DisplacementVector ) {
        type = ISVT_VECTOR;
    } else if ( valID == FluxVector ) {
        type = ISVT_SCALAR;
    } else {
        OOFEM_ERROR("unsupported UnknownType");
    }

    // print header
    if ( type == ISVT_SCALAR ) {
        fprintf(stream, "SCALARS prim_scalar_%d\n", ( int ) valID);
    } else if ( type == ISVT_VECTOR ) {
        fprintf(stream, "VECTORS vector_%d float\n", ( int ) valID);
    } else {
        OOFEM_ERROR("unsupported variable type");
    }


    SpatialLocalizer *sl = d->giveSpatialLocalizer();
    // loop over POIs
    for ( auto &poi: POIList ) {
        coords.at(1) = poi.x;
        coords.at(3) = poi.z;
        //region = poi.region;

        Element *source = sl->giveElementClosestToPoint(lcoords, closest, coords);
        if ( source ) {
            // ask interface
            EIPrimaryUnknownMapperInterface *interface =
                static_cast< EIPrimaryUnknownMapperInterface * >( source->giveInterface(EIPrimaryUnknownMapperInterfaceType) );
            if ( interface ) {
                interface->EIPrimaryUnknownMI_computePrimaryUnknownVectorAtLocal(VM_Total, tStep, lcoords, pv);
            } else {
                pv.clear();
                OOFEM_WARNING("element %d with no EIPrimaryUnknownMapperInterface support",
                               source->giveNumber() );
            }

            fprintf(stream, "%10d ", poi.id);
            if ( pv.giveSize() ) {
                for ( int j = 1; j <= pv.giveSize(); j++ ) {
                    fprintf( stream, " %15e ", pv.at(j) );
                }
            }

            fprintf(stream, "\n");
        } else {
            OOFEM_ERROR("no element containing POI(%e,%e,%e) found",
                        coords.at(1), coords.at(2), coords.at(3) );
        }
    }
}

int
EIPrimaryUnknownMapper :: evaluateAt(FloatArray &answer, IntArray &dofMask, ValueModeType mode,
                                     Domain *oldd, FloatArray &coords, IntArray &regList, TimeStep *tStep)
{
    Element *oelem;
    EIPrimaryUnknownMapperInterface *interface;
    SpatialLocalizer *sl = oldd->giveSpatialLocalizer();

    FloatArray lcoords, closest;
    if ( regList.isEmpty() ) {
        oelem = sl->giveElementClosestToPoint(lcoords, closest, coords, 0);
    } else {
        // Take the minimum of any region
        double mindist = 0.0, distance;
        oelem = NULL;
        for ( int i = 1; i < regList.giveSize(); ++i ) {
            Element *tmpelem = sl->giveElementClosestToPoint( lcoords, closest, coords, regList.at(i) );
            distance = closest.distance_square(coords);
            if ( tmpelem != NULL ) {
                distance = closest.distance_square(coords);
                if ( distance < mindist || i == 1 ) {
                    mindist = distance;
                    oelem = tmpelem;
                    if ( distance == 0.0 ) {
                        break;
                    }
                }
            }
        }
    }
    if ( !oelem ) {
        OOFEM_WARNING("Couldn't find any element containing point.");
        return false;
    }

    interface = static_cast< EIPrimaryUnknownMapperInterface * >( oelem->giveInterface(EIPrimaryUnknownMapperInterfaceType) );
    if ( interface ) {
        oelem->giveElementDofIDMask(dofMask);
        interface->EIPrimaryUnknownMI_computePrimaryUnknownVectorAtLocal(mode, tStep, lcoords, answer);
    } else {
        OOFEM_ERROR("Element does not support EIPrimaryUnknownMapperInterface");
    }

    return true;
}

void PrescribedGradientBCWeak :: assembleTangentGPContributionNew(FloatMatrix &oTangent, TracSegArray &iEl, GaussPoint &iGP, const double &iScaleFactor, const FloatArray &iBndCoord)
{
    int dim = domain->giveNumberOfSpatialDimensions();
    double detJ = 0.5 * iEl.giveLength();

    //////////////////////////////////
    // Compute traction N-matrix
    // For now, assume piecewise constant approx
    FloatArray Ntrac = FloatArray { 1.0 };
    FloatMatrix NtracMat;
    NtracMat.beNMatrixOf(Ntrac, dim);


    //////////////////////////////////
    // Compute displacement N-matrix
    // Identify the displacement element
    // we are currently standing in
    // and compute local coordinates on
    // the displacement element
    SpatialLocalizer *localizer = domain->giveSpatialLocalizer();
    FloatArray dispElLocCoord, closestPoint;
    Element *dispEl = localizer->giveElementClosestToPoint(dispElLocCoord, closestPoint, iBndCoord );


    // Compute basis functions
    XfemElementInterface *xfemElInt = dynamic_cast< XfemElementInterface * >( dispEl );
    FloatMatrix NdispMat;

    if ( xfemElInt != NULL && domain->hasXfemManager() ) {
        // If the element is an XFEM element, we use the XfemElementInterface to compute the N-matrix
        // of the enriched element.
        xfemElInt->XfemElementInterface_createEnrNmatrixAt(NdispMat, dispElLocCoord, * dispEl, false);
    } else   {
        // Otherwise, use the usual N-matrix.
        const int numNodes = dispEl->giveNumberOfDofManagers();
        FloatArray N(numNodes);

        const int dim = dispEl->giveSpatialDimension();

        NdispMat.resize(dim, dim * numNodes);
        NdispMat.zero();
        dispEl->giveInterpolation()->evalN( N, dispElLocCoord, FEIElementGeometryWrapper(dispEl) );

        NdispMat.beNMatrixOf(N, dim);
    }

    FloatMatrix contrib;
    contrib.beTProductOf(NtracMat, NdispMat);
    contrib.times( iScaleFactor * detJ * iGP.giveWeight() );

    oTangent = contrib;
}

void
POIExportModule :: exportPrimVarAs(UnknownType valID, FILE *stream, TimeStep *tStep)
{
    Domain *d = emodel->giveDomain(1);
    FloatArray pv, coords(3), lcoords, closest;
    InternalStateValueType type = ISVT_UNDEFINED;

    if ( valID == DisplacementVector ) {
        type = ISVT_VECTOR;
    } else if ( valID == FluxVector || valID == Humidity ) {
        type = ISVT_SCALAR;
    } else {
        OOFEM_ERROR("unsupported UnknownType");
    }

    // print header
    if ( type == ISVT_SCALAR ) {
        fprintf(stream, "SCALARS prim_scalar_%d\n", ( int ) valID);
    } else if ( type == ISVT_VECTOR ) {
        fprintf(stream, "VECTORS vector_%d float\n", ( int ) valID);
    } else {
        OOFEM_ERROR("unsupported variable type");
    }


    SpatialLocalizer *sl = d->giveSpatialLocalizer();
    // loop over POIs
    for ( auto &poi: POIList ) {
        coords.at(1) = poi.x;
        coords.at(3) = poi.z;
        //region = poi.region;

        Element *source = sl->giveElementClosestToPoint(lcoords, closest, coords);
        if ( source ) {
            // ask interface
            source->computeField(VM_Total, tStep, lcoords, pv);

            fprintf(stream, "%10d ", poi.id);
            for ( auto &p : pv ) {
                fprintf( stream, " %15e ", p );
            }

            fprintf(stream, "\n");
        } else {
            OOFEM_ERROR("no element containing POI(%e,%e,%e) found",
                        coords.at(1), coords.at(2), coords.at(3) );
        }
    }
}

void PrescribedGradientBCWeak :: computeIntForceGPContrib(FloatArray &oContrib_disp, IntArray &oDisp_loc_array, FloatArray &oContrib_trac, IntArray &oTrac_loc_array,TracSegArray &iEl, GaussPoint &iGP, int iDim, TimeStep *tStep, const FloatArray &iBndCoord, const double &iScaleFac, ValueModeType mode, CharType type, const UnknownNumberingScheme &s)
{

    SpatialLocalizer *localizer = domain->giveSpatialLocalizer();

	FloatMatrix contrib;
	assembleTangentGPContributionNew(contrib, iEl, iGP, iScaleFac, iBndCoord);

    // Compute vector of traction unknowns
    FloatArray tracUnknowns;
    iEl.mFirstNode->giveUnknownVector(tracUnknowns, giveTracDofIDs(), mode, tStep);

    iEl.giveTractionLocationArray(oTrac_loc_array, type, s);

    FloatArray dispElLocCoord, closestPoint;
    Element *dispEl = localizer->giveElementClosestToPoint(dispElLocCoord, closestPoint, iBndCoord );

    // Compute vector of displacement unknowns
    FloatArray dispUnknowns;
    int numDMan = dispEl->giveNumberOfDofManagers();
    for(int i = 1; i <= numDMan; i++) {
    	FloatArray nodeUnknowns;
    	DofManager *dMan = dispEl->giveDofManager(i);

    	IntArray dispIDs = giveRegularDispDofIDs();
        if(domain->hasXfemManager()) {
        	XfemManager *xMan = domain->giveXfemManager();
        	dispIDs.followedBy(xMan->giveEnrichedDofIDs(*dMan));
        }

        dMan->giveUnknownVector(nodeUnknowns, dispIDs,mode, tStep);
        dispUnknowns.append(nodeUnknowns);

    }

    dispEl->giveLocationArray(oDisp_loc_array, s);


    oContrib_disp.beTProductOf(contrib, tracUnknowns);
    oContrib_disp.negated();

    oContrib_trac.beProductOf(contrib, dispUnknowns);
    oContrib_trac.negated();
}

int
EIPrimaryUnknownMapper :: evaluateAt(FloatArray &answer, IntArray &dofMask, ValueModeType mode,
                                     Domain *oldd, FloatArray &coords, IntArray &regList, TimeStep *tStep)
{
    Element *oelem;
    EIPrimaryUnknownMapperInterface *interface;
    SpatialLocalizer *sl = oldd->giveSpatialLocalizer();

    ///@todo Change to the other version after checking that it works properly. Will render "giveElementCloseToPoint" obsolete (superseeded by giveElementClosestToPoint).
#if 1
    if ( regList.isEmpty() ) {
        oelem = sl->giveElementContainingPoint(coords);
    } else {
        oelem = sl->giveElementContainingPoint(coords, & regList);
    }
    if ( !oelem ) {
        if ( regList.isEmpty() ) {
            oelem = oldd->giveSpatialLocalizer()->giveElementCloseToPoint(coords);
        } else {
            oelem = oldd->giveSpatialLocalizer()->giveElementCloseToPoint(coords, & regList);
        }
        if ( !oelem ) {
            OOFEM_WARNING("Couldn't find any element containing point.");
            return false;
        }
    }
#else
    FloatArray lcoords, closest;
    if ( regList.isEmpty() ) {
        oelem = sl->giveElementClosestToPoint(lcoords, closest, coords, 0);
    } else {
        // Take the minimum of any region
        double mindist = 0.0, distance;
        oelem = NULL;
        for ( int i = 1; i < regList.giveSize(); ++i ) {
            Element *tmpelem = sl->giveElementClosestToPoint( lcoords, closest, coords, regList.at(i) );
            distance = closest.distance_square(coords);
            if ( tmpelem != NULL ) {
                distance = closest.distance_square(coords);
                if ( distance < mindist || i == 1 ) {
                    mindist = distance;
                    oelem = tmpelem;
                    if ( distance == 0.0 ) {
                        break;
                    }
                }
            }
        }
    }
    if ( !oelem ) {
        OOFEM_WARNING("Couldn't find any element containing point.");
        return false;
    }
#endif

    interface = static_cast< EIPrimaryUnknownMapperInterface * >( oelem->giveInterface(EIPrimaryUnknownMapperInterfaceType) );
    if ( interface ) {
        oelem->giveElementDofIDMask(dofMask);
#if 1
        FloatArray lcoords;
        if ( oelem->computeLocalCoordinates(lcoords, coords) ) {
            interface->EIPrimaryUnknownMI_computePrimaryUnknownVectorAtLocal(mode, tStep, lcoords, answer);
        } else {
            answer.clear();
        }
#else
        interface->EIPrimaryUnknownMI_computePrimaryUnknownVectorAtLocal(mode, tStep, lcoords, answer);
#endif
    } else {
        OOFEM_ERROR("Element does not support EIPrimaryUnknownMapperInterface");
    }

    return true;
}

void PrescribedGradientBCWeak :: assembleGPContrib(SparseMtrx &answer, TimeStep *tStep,
                      CharType type, const UnknownNumberingScheme &r_s, const UnknownNumberingScheme &c_s, TracSegArray &iEl, GaussPoint &iGP)
{

    SpatialLocalizer *localizer = domain->giveSpatialLocalizer();

    ///////////////
    // Gamma_plus
	FloatMatrix contrib;
	assembleTangentGPContributionNew(contrib, iEl, iGP, -1.0, iGP.giveGlobalCoordinates());

    // Compute vector of traction unknowns
    FloatArray tracUnknowns;
    iEl.mFirstNode->giveUnknownVector(tracUnknowns, giveTracDofIDs(), VM_Total, tStep);

    IntArray trac_rows;
    iEl.giveTractionLocationArray(trac_rows, type, r_s);


    FloatArray dispElLocCoord, closestPoint;
    Element *dispEl = localizer->giveElementClosestToPoint(dispElLocCoord, closestPoint, iGP.giveGlobalCoordinates() );

    IntArray disp_cols;
    dispEl->giveLocationArray(disp_cols, c_s);

    answer.assemble(trac_rows, disp_cols, contrib);

    FloatMatrix contribT;
    contribT.beTranspositionOf(contrib);
    answer.assemble(disp_cols, trac_rows, contribT);



    ///////////////
    // Gamma_minus
    contrib.clear();
	FloatArray xMinus;
	this->giveMirroredPointOnGammaMinus(xMinus, iGP.giveGlobalCoordinates() );
	assembleTangentGPContributionNew(contrib, iEl, iGP, 1.0, xMinus);

    // Compute vector of traction unknowns
	tracUnknowns.clear();
	iEl.mFirstNode->giveUnknownVector(tracUnknowns, giveTracDofIDs(), VM_Total, tStep);

    trac_rows.clear();
    iEl.giveTractionLocationArray(trac_rows, type, r_s);


    dispElLocCoord.clear(); closestPoint.clear();
    dispEl = localizer->giveElementClosestToPoint(dispElLocCoord, closestPoint, xMinus );

    disp_cols.clear();
    dispEl->giveLocationArray(disp_cols, c_s);

    answer.assemble(trac_rows, disp_cols, contrib);

    contribT.clear();
    contribT.beTranspositionOf(contrib);
    answer.assemble(disp_cols, trac_rows, contribT);


    // Assemble zeros on diagonal (required by PETSc solver)
	FloatMatrix KZero(1,1);
	KZero.zero();
    for( int i :  trac_rows) {
        answer.assemble(IntArray({i}), IntArray({i}), KZero);
    }
}

bool PLMaterialForce :: propagateInterface(Domain &iDomain, EnrichmentFront &iEnrFront, TipPropagation &oTipProp)
{
//    printf("Entering PLMaterialForce :: propagateInterface().\n");

    if ( !iEnrFront.propagationIsAllowed() ) {
        return false;
    }

    // Fetch crack tip data
    const TipInfo &tipInfo = iEnrFront.giveTipInfo();

    // Check if the tip is located in the domain
    SpatialLocalizer *localizer = iDomain.giveSpatialLocalizer();
    FloatArray lCoords, closest;
//    printf("tipInfo.mGlobalCoord: \n"); tipInfo.mGlobalCoord.printYourself();
    if( tipInfo.mGlobalCoord.giveSize() == 0 ) {
    	return false;
    }

    localizer->giveElementClosestToPoint(lCoords, closest, tipInfo.mGlobalCoord);

    if(closest.distance(tipInfo.mGlobalCoord) > 1.0e-9) {
//        printf("Tip is outside all elements.\n");
        return false;
    }


    FloatArray matForce;
    TimeStep *tStep = iDomain.giveEngngModel()->giveCurrentStep();
    mpMaterialForceEvaluator->computeMaterialForce(matForce, iDomain, tipInfo, tStep, mRadius);

//    printf("matForce: "); matForce.printYourself();

    if(matForce.giveSize() == 0) {
        return false;
    }

    double forceNorm = matForce.computeNorm();
//    printf("forceNorm: %e mCrackPropThreshold: %e\n", forceNorm, mCrackPropThreshold);

    if(forceNorm < mCrackPropThreshold || forceNorm < 1.0e-20) {
        return false;
    }

    printf("forceNorm: %e mCrackPropThreshold: %e\n", forceNorm, mCrackPropThreshold);
    printf("Propagating crack in PLMaterialForce :: propagateInterface.\n");
//    printf("Tip coord: "); tipInfo.mGlobalCoord.printYourself();

    FloatArray dir(matForce);
    dir.times(1.0/forceNorm);
//    printf("dir: "); dir.printYourself();

    const double cosAngTol = 1.0/sqrt(2.0);
    if(tipInfo.mTangDir.dotProduct(dir) < cosAngTol) {
        // Do not allow sharper turns than 45 degrees

        if( tipInfo.mNormalDir.dotProduct(dir) > 0.0 ) {
            dir = tipInfo.mTangDir;
            dir.add(tipInfo.mNormalDir);
            dir.normalize();
        }
        else {
//            dir = tipInfo.mNormalDir;
//            dir.times(-1.0);
            dir = tipInfo.mTangDir;
            dir.add(-1.0,tipInfo.mNormalDir);
            dir.normalize();
        }

        printf("//////////////////////////////////////////// Resticting crack propagation direction.\n");
//        printf("tipInfo.mTangDir: "); tipInfo.mTangDir.printYourself();
//        printf("dir: "); dir.printYourself();
    }

    // Fill up struct
    oTipProp.mTipIndex = tipInfo.mTipIndex;
    oTipProp.mPropagationDir = dir;
    oTipProp.mPropagationLength = mIncrementLength;


    return true;
}

void HangingNode :: postInitialize()
{
    Node :: postInitialize();

    Element *e;
    FEInterpolation *fei;
    FloatArray lcoords, masterContribution;

#ifdef __OOFEG
    if ( initialized ) {
        return;
    }
    initialized = true;
#endif

    // First check element and interpolation
    if ( masterElement == -1 ) { // Then we find it by taking the closest (probably containing element)
        FloatArray closest;
        SpatialLocalizer *sp = this->domain->giveSpatialLocalizer();
        sp->init();
        // Closest point or containing point? It should be contained, but with numerical errors it might be slightly outside
        // so the closest point is more robust.
        if ( !( e = sp->giveElementClosestToPoint(lcoords, closest, coordinates, this->masterRegion) ) ) {
            OOFEM_ERROR("Couldn't find closest element (automatically).");
        }
        this->masterElement = e->giveNumber();
    } else if ( !( e = this->giveDomain()->giveElement(this->masterElement) ) ) {
        OOFEM_ERROR("Requested element %d doesn't exist.", this->masterElement);
    }
    if ( !( fei = e->giveInterpolation() ) ) {
        OOFEM_ERROR("Requested element %d doesn't have a interpolator.", this->masterElement);
    }

    if ( lcoords.giveSize() == 0 ) { // we don't need to do this again if the spatial localizer was used.
        fei->global2local( lcoords, coordinates, FEIElementGeometryWrapper(e) );
    }

    // Initialize slave dofs (inside check of consistency of receiver and master dof)
    const IntArray &masterNodes = e->giveDofManArray();
    for ( Dof *dof: *this ) {
        SlaveDof *sdof = dynamic_cast< SlaveDof * >(dof);
        if ( sdof ) {
            DofIDItem id = sdof->giveDofID();
            fei = e->giveInterpolation(id);
            if ( !fei ) {
                OOFEM_ERROR("Requested interpolation for dof id %d doesn't exist in element %d.",
                             id, this->masterElement);
            }
#if 0 // This won't work (yet), as it requires some more general FEI classes, or something similar.
            if ( fei->hasMultiField() ) {
                FloatMatrix multiContribution;
                IntArray masterDofIDs, masterNodesDup, dofids;
                fei->evalMultiN(multiContribution, dofids, lcoords, FEIElementGeometryWrapper(e), 0.0);
                masterContribution.flatten(multiContribution);
                masterDofIDs.clear();
                for ( int i = 0; i <= multiContribution.giveNumberOfColumns(); ++i ) {
                    masterDofIDs.followedBy(dofids);
                    masterNodesDup.followedBy(masterNodes);
                }
                sdof->initialize(masterNodesDup, & masterDofIDs, masterContribution);
            } else { }
#else
            // Note: There can be more masterNodes than masterContributions, since all the
            // FEI classes are based on that the first nodes correspond to the simpler/linear interpolation.
            // If this assumption is changed in FEIElementGeometryWrapper + friends,
            // masterNode will also need to be modified for each dof accordingly.
            fei->evalN( masterContribution, lcoords, FEIElementGeometryWrapper(e) );
            sdof->initialize(masterNodes, IntArray(), masterContribution);
#endif
        }
    }
}