C++ SpatialLocalizer类代码示例

void
MMAContainingElementProjection :: __init(Domain *dold, IntArray &type, FloatArray &coords, Set &elemSet, TimeStep *tStep, bool iCohesiveZoneGP)
{
    SpatialLocalizer *sl = dold->giveSpatialLocalizer();
    FloatArray jGpCoords;
    double distance, minDist = 1.e6;
    Element *srcElem;

    if ( ( srcElem = sl->giveElementContainingPoint(coords, elemSet) ) ) {
        this->source = NULL;
        for ( GaussPoint *jGp: *srcElem->giveDefaultIntegrationRulePtr() ) {
            if ( srcElem->computeGlobalCoordinates( jGpCoords, jGp->giveNaturalCoordinates() ) ) {
                distance = coords.distance(jGpCoords);
                if ( distance < minDist ) {
                    minDist = distance;
                    this->source = jGp;
                }
            }
        }

        if ( !source ) {
            OOFEM_ERROR("no suitable source found");
        }
    } else {
        OOFEM_ERROR("No suitable element found");
    }
}

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

bool PLCrackPrescribedDir :: propagateInterface(Domain &iDomain, EnrichmentFront &iEnrFront, TipPropagation &oTipProp)
{
    if ( !iEnrFront.propagationIsAllowed() ) {
        return false;
    }

    const TipInfo &tipInfo = iEnrFront.giveTipInfo();

    SpatialLocalizer *localizer = iDomain.giveSpatialLocalizer();
    // It is meaningless to propagate a tip that is not inside any element
    if ( tipInfo.mGlobalCoord.giveSize() == 0 ) {
        return false;
    }

    Element *el = localizer->giveElementContainingPoint(tipInfo.mGlobalCoord);
    if ( el == NULL ) {
        return false;
    }


    double angleRad = mAngle * M_PI / 180.0;
    FloatArray dir = {
        cos(angleRad), sin(angleRad)
    };

    oTipProp.mTipIndex = tipInfo.mTipIndex;
    oTipProp.mPropagationDir = dir;
    oTipProp.mPropagationLength = mIncrementLength;

    return true;
}

int
PrimaryField :: __evaluateAt(FloatArray &answer, FloatArray& coords,
                    ValueModeType mode, TimeStep *atTime,
                    IntArray *dofId)
{
    Element *bgelem;
    Domain *domain = emodel->giveDomain(domainIndx);
    SpatialLocalizer *sl = domain->giveSpatialLocalizer();
    // locate background element
    if ( ( bgelem = sl->giveElementContainingPoint(coords) ) == NULL ) {
        //_error ("PrimaryField::evaluateAt: point not found in domain\n");
        return 1;
    }

    EIPrimaryFieldInterface *interface = ( EIPrimaryFieldInterface * ) ( bgelem->giveInterface(EIPrimaryFieldInterfaceType) );
    if ( interface ) {
        if (dofId) {
            return interface->EIPrimaryFieldI_evaluateFieldVectorAt(answer, * this, coords, *dofId, mode, atTime);
        } else { // use element default dof id mask
            IntArray elemDofId;
            bgelem->giveElementDofIDMask(this->giveEquationID(), elemDofId);
            return interface->EIPrimaryFieldI_evaluateFieldVectorAt(answer, * this, coords, elemDofId, mode, atTime);
        }
    } else {
        _error("ScalarPrimaryField::operator(): background element does not support EIPrimaryFiledInterface\n");
        return 1; // failed
    }
}

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

void
Delamination :: findInitiationFronts(bool &failureChecked, const IntArray &CSnumbers, std :: vector< IntArray > &CSinterfaceNumbers, std :: vector< IntArray > &CSDofManNumbers, std :: vector< FloatArray > &initiationFactors, TimeStep *tStep)
{           
    // Loop through all cross sections associated with delaminations 
    // Returns 
    // CSinterfaceNumbers: the failed interface number associated with the cross sections. 
    // CSDofManNumbers: the dofmanagers that should be enriched (associated with the cross sections)
    IntArray failedElementInterfaces;  
    IntArray elementNumbers; 
    SpatialLocalizer *localizer = this->giveDomain()->giveSpatialLocalizer();
    
    // NB: Assumes that elements can only be included in one cross section. 
    
    for ( int iCS = 1 ; iCS <= CSnumbers.giveSize() ; iCS++ ) {
        
        int eltSetNumber = this->giveDomain()->giveCrossSection(CSnumbers.at(iCS))->giveSetNumber();
        //printf("Cross section No. %i, set No. %i \n",CSnumbers.at(iCS),eltSetNumber);
        IntArray elementNumbers = this->giveDomain()->giveSet(eltSetNumber)->giveElementList();
        
        for ( auto eltNumber : elementNumbers ) {
            Element *elt = this->giveDomain()->giveGlobalElement(eltNumber);
            
            if ( Shell7BaseXFEM *shellElt = dynamic_cast < Shell7BaseXFEM * > (elt) ) {
                
                //bool recoverStresses = true;
                shellElt->giveFailedInterfaceNumber(failedElementInterfaces, initiationFactors[iCS-1], tStep, this->recoverStresses);
                //failedElementInterfaces.printYourself("failedElementInterfaces");
                for (int eltInt : failedElementInterfaces ) {
                    CSinterfaceNumbers[iCS-1].insertSortedOnce(eltInt);
                }
                if ( !failedElementInterfaces.isEmpty() ) {
                    for (int iDF : shellElt->giveDofManArray() ) {
                        //printf("element node %d \n",iDF);
                        if ( this->initiationRadius > 0.0 ) {
                            const FloatArray gCoords = this->giveDomain()->giveNode(iDF)->giveNodeCoordinates();
                            std :: list< int > nodeList;
                            localizer->giveAllNodesWithinBox(nodeList,gCoords,initiationRadius);
                                
                            for ( int jNode : nodeList ) {
                                //printf("nodeList node %d \n",jNode);
                                CSDofManNumbers[iCS-1].insertSortedOnce(jNode);
                            }
                        } else {
                            CSDofManNumbers[iCS-1].insertSortedOnce(iDF);
                        }
                    }
                }
            }
        }
    } 
    
    failureChecked = 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
MMAClosestIPTransfer :: __init(Domain *dold, IntArray &type, FloatArray &coords, Set &elemSet, TimeStep *tStep, bool iCohesiveZoneGP)
{
    SpatialLocalizer *sl = dold->giveSpatialLocalizer();
    this->source = sl->giveClosestIP(coords, elemSet, iCohesiveZoneGP);

    if ( !source ) {
        OOFEM_ERROR("no suitable source found");
    }

    mpMaterialStatus = dynamic_cast<MaterialStatus*>(source->giveMaterialStatus());
    if( mpMaterialStatus == NULL ) {
        OOFEM_ERROR("Could not find material status.");
    }
}

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

void PrescribedGradientBCPeriodic :: findSlaveToMasterMap()
{
    FloatArray coord;
    SpatialLocalizer *sl = this->domain->giveSpatialLocalizer();
    //Set *masterSet = this->domain->giveSet(2);
    const IntArray &nodes = this->domain->giveSet(this->set)->giveNodeList(); // Split into slave set and master set?
    int nsd = jump.giveSize();
    std :: vector< FloatArray > jumps;
    // Construct all the possible jumps;
    jumps.reserve((2 << (nsd-1)) - 1);
    if ( nsd != 3 ) {
        OOFEM_ERROR("Only 3d implemented yet!");
    }
    jumps.emplace_back(jump);
    jumps.emplace_back(FloatArray{jump.at(1), jump.at(2), 0.});
    jumps.emplace_back(FloatArray{jump.at(1), 0., jump.at(3)});
    jumps.emplace_back(FloatArray{0., jump.at(2), jump.at(3)});
    jumps.emplace_back(FloatArray{jump.at(1), 0., 0.});
    jumps.emplace_back(FloatArray{0., jump.at(2), 0.});
    jumps.emplace_back(FloatArray{0., 0., jump.at(3)});

    this->slavemap.clear();
    for ( int inode : nodes ) {
        Node *masterNode = NULL;
        Node *node = this->domain->giveNode(inode);
        const FloatArray &masterCoord = *node->giveCoordinates();
        //printf("node %d\n", node->giveLabel()); masterCoord.printYourself();
        // The difficult part, what offset to subtract to find the master side;
        for ( FloatArray &testJump : jumps ) {
            coord.beDifferenceOf(masterCoord, testJump);
            masterNode = sl->giveNodeClosestToPoint(coord, fabs(jump.at(1))*1e-5);
            if ( masterNode != NULL ) {
                //printf("Found master (%d) to node %d (distance = %e)\n",  masterNode->giveNumber(), node->giveNumber(),
                //       masterNode->giveCoordinates()->distance(coord));
                break;
            }
        }
        if ( masterNode != NULL ) {
            this->slavemap.insert({node->giveNumber(), masterNode->giveNumber()});
        } else {
            OOFEM_ERROR("Couldn't find master node!");
        }
    }
}

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
MMALeastSquareProjection :: __init(Domain *dold, IntArray &type, FloatArray &coords, Set &elemSet, TimeStep *tStep, bool iCohesiveZoneGP)
//(Domain* dold, IntArray& varTypes, GaussPoint* gp, TimeStep* tStep)
{
    GaussPoint *sourceIp;
    Element *sourceElement;
    SpatialLocalizer *sl = dold->giveSpatialLocalizer();
    IntegrationRule *iRule;

    IntArray patchList;

    this->patchDomain = dold;
    // find the closest IP on old mesh
    sourceElement = sl->giveElementContainingPoint(coords, elemSet);

    if ( !sourceElement ) {
        OOFEM_ERROR("no suitable source element found");
    }

    // determine the type of patch
    Element_Geometry_Type egt = sourceElement->giveGeometryType();
    if ( egt == EGT_line_1 ) {
        this->patchType = MMALSPPatchType_1dq;
    } else if ( ( egt == EGT_triangle_1 ) || ( egt == EGT_quad_1 ) ) {
        this->patchType = MMALSPPatchType_2dq;
    } else {
        OOFEM_ERROR("unsupported material mode");
    }

    /* Determine the state of closest point.
     * Only IP in the neighbourhood with same state can be used
     * to interpolate the values.
     */
    FloatArray dam;
    int state = 0;
    if ( this->stateFilter ) {
        iRule = sourceElement->giveDefaultIntegrationRulePtr();
        for ( GaussPoint *gp: *iRule ) {
            sourceElement->giveIPValue(dam, gp, IST_PrincipalDamageTensor, tStep);
            if ( dam.computeNorm() > 1.e-3 ) {
                state = 1; // damaged
            }
        }
    }

    // from source neighbours the patch will be constructed
    Element *element;
    IntArray neighborList;
    patchList.resize(1);
    patchList.at(1) = sourceElement->giveNumber();
    int minNumberOfPoints = this->giveNumberOfUnknownPolynomialCoefficients(this->patchType);
    int actualNumberOfPoints = sourceElement->giveDefaultIntegrationRulePtr()->giveNumberOfIntegrationPoints();
    int nite = 0;
    int elemFlag;
    // check if number of IP in patchList is sufficient
    // some recursion control would be appropriate
    while ( ( actualNumberOfPoints < minNumberOfPoints ) && ( nite <= 2 ) ) {
        //if not,  construct the neighborhood
        dold->giveConnectivityTable()->giveElementNeighbourList(neighborList, patchList);
        // count number of available points
        patchList.clear();
        actualNumberOfPoints = 0;
        for ( int i = 1; i <= neighborList.giveSize(); i++ ) {
            if ( this->stateFilter ) {
                element = patchDomain->giveElement( neighborList.at(i) );
                // exclude elements in different regions
                if ( !elemSet.hasElement( element->giveNumber() ) ) {
                    continue;
                }

                iRule = element->giveDefaultIntegrationRulePtr();
                elemFlag = 0;
                for ( GaussPoint *gp: *iRule ) {
                    element->giveIPValue(dam, gp, IST_PrincipalDamageTensor, tStep);
                    if ( state && ( dam.computeNorm() > 1.e-3 ) ) {
                        actualNumberOfPoints++;
                        elemFlag = 1;
                    } else if ( ( state == 0 ) && ( dam.computeNorm() < 1.e-3 ) ) {
                        actualNumberOfPoints++;
                        elemFlag = 1;
                    }
                }

                if ( elemFlag ) {
                    // include this element with corresponding state in neighbor search.
                    patchList.followedBy(neighborList.at(i), 10);
                }
            } else { // if (! yhis->stateFilter)
                element = patchDomain->giveElement( neighborList.at(i) );
                // exclude elements in different regions
                if ( !elemSet.hasElement( element->giveNumber() ) ) {
                    continue;
                }

                actualNumberOfPoints += element->giveDefaultIntegrationRulePtr()->giveNumberOfIntegrationPoints();

                patchList.followedBy(neighborList.at(i), 10);
            }
        } // end loop over neighbor list

//.........这里部分代码省略.........