plprincipalstrain.C

Go to the documentation of this file.

/*
 *
 *                 #####    #####   ######  ######  ###   ###
 *               ##   ##  ##   ##  ##      ##      ## ### ##
 *              ##   ##  ##   ##  ####    ####    ##  #  ##
 *             ##   ##  ##   ##  ##      ##      ##     ##
 *            ##   ##  ##   ##  ##      ##      ##     ##
 *            #####    #####   ##      ######  ##     ##
 *
 *
 *             OOFEM : Object Oriented Finite Element Code
 *
 *               Copyright (C) 1993 - 2013   Borek Patzak
 *
 *
 *
 *       Czech Technical University, Faculty of Civil Engineering,
 *   Department of Structural Mechanics, 166 29 Prague, Czech Republic
 *
 *  This library is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU Lesser General Public
 *  License as published by the Free Software Foundation; either
 *  version 2.1 of the License, or (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *  Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public
 *  License along with this library; if not, write to the Free Software
 *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
 */

#include "plprincipalstrain.h"
#include "xfem/propagationlaw.h"
#include "xfem/tipinfo.h"
#include "classfactory.h"
#include "mathfem.h"
#include "dynamicinputrecord.h"
#include "spatiallocalizer.h"
#include "floatmatrix.h"
#include "gausspoint.h"
#include "xfem/enrichmentitem.h"
#include "feinterpol.h"
#include "xfem/xfemmanager.h"

#include "Materials/structuralms.h"
#include "Materials/structuralmaterial.h"

#include "xfem/XFEMDebugTools.h"

namespace oofem {
REGISTER_PropagationLaw(PLPrincipalStrain)

PLPrincipalStrain::PLPrincipalStrain():
mRadius(0.0),mIncrementLength(0.0),mStrainThreshold(0.0), mUseRadialBasisFunc(false)
{

}

PLPrincipalStrain::~PLPrincipalStrain() {

}

IRResultType PLPrincipalStrain :: initializeFrom(InputRecord *ir)
{
    IRResultType result;

    IR_GIVE_FIELD(ir, mRadius,                          _IFT_PLPrincipalStrain_Radius);
    IR_GIVE_FIELD(ir, mIncrementLength,         _IFT_PLPrincipalStrain_IncLength);
    IR_GIVE_FIELD(ir, mStrainThreshold, _IFT_PLPrincipalStrain_StrainThreshold);

    int useRadialBasisFunc = 0;
    IR_GIVE_OPTIONAL_FIELD(ir, useRadialBasisFunc, _IFT_PLPrincipalStrain_RadialBasisFunc);
    if ( useRadialBasisFunc == 1 ) {
        mUseRadialBasisFunc = true;
    }

    return IRRT_OK;
}

void PLPrincipalStrain :: giveInputRecord(DynamicInputRecord &input)
{
    int number = 1;
    input.setRecordKeywordField(this->giveInputRecordName(), number);

    input.setField(mRadius,                             _IFT_PLPrincipalStrain_Radius);
    input.setField(mIncrementLength,            _IFT_PLPrincipalStrain_IncLength);
    input.setField(mStrainThreshold,        _IFT_PLPrincipalStrain_StrainThreshold);

    if ( mUseRadialBasisFunc ) {
        input.setField(1,       _IFT_PLPrincipalStrain_RadialBasisFunc);
    }
}


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

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

    SpatialLocalizer *localizer = iDomain.giveSpatialLocalizer();


    const FloatArray &xT    = tipInfo.mGlobalCoord;
    const FloatArray &t     = tipInfo.mTangDir;
//    const FloatArray &n     = tipInfo.mNormalDir;

    // It is meaningless to propagate a tip that is not inside any element
    Element *el = localizer->giveElementContainingPoint(tipInfo.mGlobalCoord);
    if ( el != NULL ) {



      FloatArray x(xT);
      x.add(mRadius, t);
//      circPoints.push_back(x);


        std :: vector< double >sigTTArray, sigRTArray;

        FloatArray strainVec;

        if ( mUseRadialBasisFunc ) {
          // Interpolate strain 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;
          IntArray elIndices;
          localizer->giveAllElementsWithIpWithinBox(elIndices, x, 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);

            for ( GaussPoint *gp_i: *gpEl->giveDefaultIntegrationRulePtr() ) {
              // Compute global gp coordinates
              FloatArray N;
              FEInterpolation *interp = gpEl->giveInterpolation();
              interp->evalN( N, gp_i->giveNaturalCoordinates(), 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;
              if ( tipToGP.dotProduct(t) < 0.0 ) {
                inFrontOfCrack = false;
              }

              double r = x.distance(globalCoord);

              if ( r < l && inFrontOfCrack ) {
                double w = ( ( l - r ) / ( pow(2.0 * M_PI, 1.5) * pow(l, 3) ) ) * exp( -0.5 * pow(r, 2) / pow(l, 2) );

                // Compute gp volume
                double V = gpEl->computeVolumeAround(gp_i);

                // Get stress
                StructuralMaterialStatus *ms = dynamic_cast< StructuralMaterialStatus * >( gp_i->giveMaterialStatus() );
                if ( ms == NULL ) {
                  OOFEM_ERROR("failed to fetch MaterialStatus.");
                }

                FloatArray strainVecGP = ms->giveStrainVector();

                if ( sumQiWiVi.giveSize() != strainVecGP.giveSize() ) {
                  sumQiWiVi.resize( strainVecGP.giveSize() );
                  sumQiWiVi.zero();
                }

                // Add to numerator
                sumQiWiVi.add(w * V, strainVecGP);

                // Add to denominator
                sumWiVi += w * V;
              }
            }
          }


          if ( fabs(sumWiVi) > 1.0e-12 ) {
            strainVec.beScaled(1.0 / sumWiVi, sumQiWiVi);
          } else {
            // Take strain from closest Gauss point
            int region = 1;
            bool useCZGP = false;
            GaussPoint &gp = * ( localizer->giveClosestIP(x, region, useCZGP) );


            // Compute strain
            StructuralMaterialStatus *ms = dynamic_cast< StructuralMaterialStatus * >( gp.giveMaterialStatus() );
            if ( ms == NULL ) {
              OOFEM_ERROR("failed to fetch MaterialStatus.");
            }

            strainVec = ms->giveStrainVector();
          }
        } else {
          // Take stress from closest Gauss point
          int region = 1;
          bool useCZGP = false;
          GaussPoint &gp = * ( localizer->giveClosestIP(x, region, useCZGP) );


          // Compute stresses
          StructuralMaterialStatus *ms = dynamic_cast< StructuralMaterialStatus * >( gp.giveMaterialStatus() );
          if ( ms == NULL ) {
            OOFEM_ERROR("failed to fetch MaterialStatus.");
          }

          strainVec = ms->giveStrainVector();
        }

        // Compute principal strain
    FloatArray principalVals;
    FloatMatrix principalDirs;
    StructuralMaterial::computePrincipalValDir(principalVals, principalDirs, strainVec, principal_strain);



        // Compare with threshold
        printf("Max principal strain: %e\n", principalVals[0]);
        if ( principalVals[0] > mStrainThreshold ) {

      FloatArray propNormal;
      propNormal.beColumnOf(principalDirs, 1);

      FloatArray propTangent = {-propNormal(1), propNormal(0)};

      if( propTangent.dotProduct(t) < 0.0 ) {
        propTangent.times(-1.0);
      }

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

            return true;
        }

    } // Tip is inside an element.


    return false;
}

} // end namespace oofem