PlateModel< T >::Impl Class Reference

#include <PlateModel_impl.h>

Collaboration diagram for PlateModel< T >::Impl:

Public Types
typedef T	T_Scalar
typedef Array< T >	TArray
typedef Vector< T, 3 >	VectorT3
typedef VectorTranspose< T, 3 >	VectorT3T
typedef Array< VectorT3 >	VectorT3Array
typedef Vector< T, 4 >	VectorT4
typedef Array< VectorT4 >	VectorT4Array
typedef SquareTensor< T, 3 >	DiagTensorT3
typedef CRMatrix< DiagTensorT3, DiagTensorT3, VectorT3 >	VVMatrix
typedef VVMatrix::DiagArray	VVDiagArray
typedef Gradient< VectorT3 >	VGradType
typedef Array< Gradient < VectorT3 > >	VGradArray

Public Member Functions
	Impl (const GeomFields &geomFields, PlateFields &plateFields, const MeshList &meshes)
void	init ()
PlateBCMap &	getBCMap ()
PlateVCMap &	getVCMap ()
PlateModelOptions< T > &	getOptions ()
void	updateTime ()
void	initDeformationLinearization (LinearSystem &ls)
void	linearizeDeformation (LinearSystem &ls)
MFRPtr	solveDeformation ()
void	calculatePlateVelocity ()
void	calculatePlateAcceleration ()
void	postPlateSolve (LinearSystem &ls)
bool	advance (const int niter)
void	printBCs ()
void	getMoment (const Mesh &mesh)
map< string, shared_ptr < ArrayBase > > &	getPersistenceData ()
void	restart ()
void	dumpMatrix (const string fileBase)

Private Attributes
const MeshList	_meshes
const GeomFields &	_geomFields
PlateFields &	_plateFields
PlateBCMap	_bcMap
PlateVCMap	_vcMap
PlateModelOptions< T >	_options
GradientModel< VectorT3 >	_deformationGradientModel
MFRPtr	_initialDeformationNorm
int	_niters
map< string, shared_ptr < ArrayBase > >	_persistenceData

Detailed Description

template<class T>
class PlateModel< T >::Impl

Definition at line 298 of file PlateModel_impl.h.

Member Typedef Documentation

template<class T>

typedef SquareTensor<T,3> PlateModel< T >::Impl::DiagTensorT3

Definition at line 311 of file PlateModel_impl.h.

template<class T>

typedef T PlateModel< T >::Impl::T_Scalar

Definition at line 301 of file PlateModel_impl.h.

template<class T>

typedef Array<T> PlateModel< T >::Impl::TArray

Definition at line 302 of file PlateModel_impl.h.

template<class T>

typedef Vector<T,3> PlateModel< T >::Impl::VectorT3

Definition at line 303 of file PlateModel_impl.h.

template<class T>

typedef Array<VectorT3> PlateModel< T >::Impl::VectorT3Array

Definition at line 306 of file PlateModel_impl.h.

template<class T>

typedef VectorTranspose<T,3> PlateModel< T >::Impl::VectorT3T

Definition at line 304 of file PlateModel_impl.h.

template<class T>

typedef Vector<T,4> PlateModel< T >::Impl::VectorT4

Definition at line 307 of file PlateModel_impl.h.

template<class T>

typedef Array<VectorT4> PlateModel< T >::Impl::VectorT4Array

Definition at line 308 of file PlateModel_impl.h.

template<class T>

typedef Array<Gradient<VectorT3> > PlateModel< T >::Impl::VGradArray

Definition at line 318 of file PlateModel_impl.h.

template<class T>

typedef Gradient<VectorT3> PlateModel< T >::Impl::VGradType

Definition at line 317 of file PlateModel_impl.h.

template<class T>

typedef VVMatrix::DiagArray PlateModel< T >::Impl::VVDiagArray

Definition at line 314 of file PlateModel_impl.h.

template<class T>

typedef CRMatrix<DiagTensorT3,DiagTensorT3,VectorT3> PlateModel< T >::Impl::VVMatrix

Definition at line 313 of file PlateModel_impl.h.

Constructor & Destructor Documentation

template<class T>

PlateModel< T >::Impl::Impl ( const GeomFields &  geomFields, PlateFields &  plateFields, const MeshList &  meshes  )

inline

Definition at line 320 of file PlateModel_impl.h.

References Model::_meshes, PlateBC< T >::bcType, Mesh::getAllFaceGroups(), Mesh::getID(), FaceGroup::groupType, FaceGroup::id, and PlateVC< T >::vcType.

                               :
    _meshes(meshes),
    _geomFields(geomFields),
    _plateFields(plateFields),
    _deformationGradientModel(_meshes,_plateFields.deformation,
                           _plateFields.deformationGradient,_geomFields),
    _initialDeformationNorm(),
    _niters(0)
  {
    const int numMeshes = _meshes.size();
    for (int n=0; n<numMeshes; n++)
    {
        const Mesh& mesh = *_meshes[n];

        PlateVC<T> *vc(new PlateVC<T>());
        vc->vcType = "plate";
        _vcMap[mesh.getID()] = vc;
        foreach(const FaceGroupPtr fgPtr, mesh.getAllFaceGroups())
        {
            const FaceGroup& fg = *fgPtr;
            if ((_bcMap.find(fg.id) == _bcMap.end())&&(fg.groupType != "interior"))
            {
                PlateBC<T> *bc(new PlateBC<T>());
                
                _bcMap[fg.id] = bc;
                if (fg.groupType == "wall") 
                {
                    bc->bcType = "SpecifiedTraction";
                }
                else if (fg.groupType == "interface")
                  {
                    bc->bcType = "Interface";
                  }
                else if (fg.groupType == "symmetry")
                {
                    bc->bcType = "Symmetry";
                }
                else if ((fg.groupType == "velocity-inlet") ||
                         (fg.groupType == "pressure-inlet")) 
                {
                    bc->bcType = "SpecifiedDeformation";
                }
                else if (fg.groupType == "pressure-outlet") 
                {
                    bc->bcType = "SpecifiedForce";
                }
                else
                  throw CException("PlateModel: unknown face group type "
                                   + fg.groupType);
            }
        }
    }
  }

Member Function Documentation

template<class T>

bool PlateModel< T >::Impl::advance ( const int  niter )

inline

Definition at line 1006 of file PlateModel_impl.h.

  {

    for(int n=0; n<niter; n++)
    { 
        MFRPtr dNorm = solveDeformation();

        if (_niters < 5)
        {
            _initialDeformationNorm->setMax(*dNorm);
        }
        
        MFRPtr dNormRatio(dNorm->normalize(*_initialDeformationNorm));
        
        if (_options.printNormalizedResiduals)
          cout << _niters << ": " << *dNormRatio <<  endl;
        else
          cout << _niters << ": " << *dNorm <<  endl;

        _niters++;
        if (*dNormRatio < _options.deformationTolerance)
          return true;
    }
    return false;
  }

template<class T>

void PlateModel< T >::Impl::calculatePlateAcceleration ( )

inline

Definition at line 888 of file PlateModel_impl.h.

References Model::_meshes, Mesh::getCells(), and StorageSite::getCountLevel1().

  {

    T_Scalar two(2.0);
    T_Scalar three(3.0);
    T_Scalar five(5.0);
    T_Scalar four(4.0);
    T_Scalar twelve(12.0);
    T_Scalar dT = _options["timeStep"];
    
    const T_Scalar dT2 = dT*dT;

    const int numMeshes = _meshes.size();
    for(int n=0;n<numMeshes;n++)
    {
        const Mesh& mesh = *_meshes[n];
        const StorageSite& cells = mesh.getCells();
        const int nCells = cells.getCountLevel1();

        TArray& acceleration =
          dynamic_cast<TArray&>(_plateFields.acceleration[cells]); 
        const VectorT3Array& w =
          dynamic_cast<const VectorT3Array&>(_plateFields.deformation[cells]);
        const VectorT3Array& wN1 =
          dynamic_cast<const VectorT3Array&>(_plateFields.deformationN1[cells]);
        const VectorT3Array& wN2 =
          dynamic_cast<const VectorT3Array&>(_plateFields.deformationN2[cells]);
        const TArray& density = 
          dynamic_cast<const TArray&> (_plateFields.density[cells]);
        //constant time step
        if(!_options.variableTimeStep)
          {
            if (_options.timeDiscretizationOrder > 1)        //second order
              {
                const VectorT3Array& wN3 =
                  dynamic_cast<const VectorT3Array&>(_plateFields.deformationN3[cells]);
                for(int c=0; c<nCells; c++)
                  {                 
                    const T_Scalar rhobydT2 = density[c]/dT2;
                    acceleration[c] = rhobydT2*(two*w[c][2] - five*wN1[c][2] + four*wN2[c][2]
                                                - wN3[c][2]);
                  }
              }
            else                                              // first order
              {
                for(int c=0; c<nCells; c++)
                  {
                    const T_Scalar rhobydT2 = density[c]/dT2;
                    acceleration[c] = rhobydT2*(w[c][2]- two*wN1[c][2] + wN2[c][2]);
                  }
              }
          }
        //variable time step
        else
          {
            T_Scalar dTN1 = _options.timeStepN1;
            T_Scalar dTN2 = _options.timeStepN2;
            T_Scalar a = (dT + dTN1)/dT;
            T_Scalar b = (dT + dTN1 + dTN2)/dT;
            T_Scalar one(1.0);
            if (_options.timeDiscretizationOrder > 1)        //second order
              { 
                const VectorT3Array& wN3 =
                  dynamic_cast<const VectorT3Array&>(_plateFields.deformationN3[cells]);
                T_Scalar c1 = (two*a*b*(pow(a,two)-pow(b,two))+two*b*(pow(b,two)-one)-two*a*(pow(a,two)-one))/
                  (a*b*(a-one)*(b-one)*(a-b));
                T_Scalar c2 = -two*(a+b)/((a-1)*(b-1));
                T_Scalar c3 = -two*(b+one)/(a*(a-b)*(a-one));
                T_Scalar c4 = two*(a+one)/(b*(a-b)*(b-one));            
                for(int c=0; c<nCells; c++)
                  {
                    const T_Scalar rhobydT2 = density[c]/dT2;
                    acceleration[c] = rhobydT2*(c1*w[c][2] + c2*wN1[c][2] + c3*wN2[c][2]
                                              + c4*wN3[c][2]);
                  }
              }
            else
              {
                T_Scalar c1 = two/a;
                T_Scalar c2 = -two/(a-one);
                T_Scalar c3 = two/(a*(a-one));
                for(int c=0; c<nCells; c++)
                  {
                    const T_Scalar rhobydT2 = density[c]/dT2;
                    acceleration[c] = rhobydT2*(c1*w[c][2] + c2*wN1[c][2]
                                                + c3*wN2[c][2]);
                  }
              }
          }
    }
  }

template<class T>

void PlateModel< T >::Impl::calculatePlateVelocity ( )

inline

Definition at line 865 of file PlateModel_impl.h.

References Model::_meshes, Mesh::getCells(), and StorageSite::getCountLevel1().

  {
    const int numMeshes = _meshes.size();
    for(int n=0;n<numMeshes;n++)
    {
        const Mesh& mesh = *_meshes[n];
        const StorageSite& cells = mesh.getCells();
        const int nCells = cells.getCountLevel1();
        const VectorT3Array& w =
          dynamic_cast<const VectorT3Array&>(_plateFields.deformation[cells]);
        const VectorT3Array& wN1 =
          dynamic_cast<const VectorT3Array&>(_plateFields.deformationN1[cells]);
        VectorT3Array& velocity =
          dynamic_cast<VectorT3Array&>(_plateFields.velocity[cells]);
        const T timeStep = _options["timeStep"];
        for (int c=0; c<nCells; c++){
            velocity[c] = (w[c]-wN1[c])/timeStep;
        }
        
    }
  }

template<class T>

void PlateModel< T >::Impl::dumpMatrix ( const string  fileBase )

inline

Definition at line 1180 of file PlateModel_impl.h.

References Model::_meshes, LinearSystem::getB(), Mesh::getCells(), CRConnectivity::getCol(), CRMatrix< T_Diag, T_OffDiag, X >::getConnectivity(), CRMatrix< T_Diag, T_OffDiag, X >::getDiag(), LinearSystem::getMatrix(), MultiFieldMatrix::getMatrix(), CRMatrix< T_Diag, T_OffDiag, X >::getOffDiag(), CRConnectivity::getRow(), StorageSite::getSelfCount(), LinearSystem::initAssembly(), and LinearSystem::initSolve().

  {
    LinearSystem ls;
    initDeformationLinearization(ls);
    ls.initAssembly();
    linearizeDeformation(ls);
    ls.initSolve();
    
    MultiFieldMatrix& mfmatrix = ls.getMatrix();
    MultiField& b = ls.getB();

    const Mesh& mesh = *_meshes[0];
    const StorageSite& cells = mesh.getCells();
    
    MultiField::ArrayIndex vIndex(&_plateFields.deformation,&cells);

    VVMatrix& matrix =
      dynamic_cast<VVMatrix&>(mfmatrix.getMatrix(vIndex,vIndex));

    VVDiagArray& diag = matrix.getDiag();
    VVDiagArray& coeff = matrix.getOffDiag();

    VectorT3Array& rCell = dynamic_cast<VectorT3Array&>(b[vIndex]);

    const CRConnectivity& cr = matrix.getConnectivity();

    const Array<int>& row = cr.getRow();
    const Array<int>& col = cr.getCol();
    
    const int nCells = cells.getSelfCount();

    const int dimension = 3;

    const int dim2 = dimension*dimension;
    
    int nFlatRows = dimension*nCells;

    int nFlatCoeffs = nCells*dim2;

    for(int i=0; i<nCells; i++)
      for(int jp=row[i]; jp<row[i+1]; jp++)
      {
          const int j = col[jp];
          if (j<nCells) nFlatCoeffs += dim2;
      }
    
    string matFileName = fileBase + ".mat";
    FILE *matFile = fopen(matFileName.c_str(),"wb");
    
    fprintf(matFile,"%%%%MatrixMarket matrix coordinate real general\n");
    fprintf(matFile,"%d %d %d\n", nFlatRows,nFlatRows,nFlatCoeffs);

    for(int i=0; i<nCells; i++)
    {
        for(int ndr=0; ndr<dimension; ndr++)
          for(int ndc=0; ndc<dimension; ndc++)
          {
              const int nfr = i*dimension + ndr;
              const int nfc = i*dimension + ndc;
              T ap = diag[i](ndr,ndc);
              //if (fabs(ap) < 1.0) ap =0.;
              fprintf(matFile,"%d %d %22.15le\n", nfr+1, nfc+1, ap);
          }
        
        for(int jp=row[i]; jp<row[i+1]; jp++)
        {
            const int j = col[jp];
            if (j<nCells)
            {
                for(int ndr=0; ndr<dimension; ndr++)
                  for(int ndc=0; ndc<dimension; ndc++)
                  {
                      const int nfr = i*dimension + ndr;
                      const int nfc = j*dimension + ndc;
                      T anb = coeff[jp](ndr,ndc);
                      //if (fabs(anb) < 1.0) anb =0.;
                      fprintf(matFile,"%d %d %22.15le\n", nfr+1, nfc+1, anb);
                  }
            }
        }
    }
    
    fclose(matFile);

    string rhsFileName = fileBase + ".rhs";
    FILE *rhsFile = fopen(rhsFileName.c_str(),"wb");
    
    for(int i=0; i<nCells; i++)
      for(int nd=0; nd < dimension; nd++)
            fprintf(rhsFile,"%22.15le\n",-rCell[i][nd]);

    fclose(rhsFile);
  }

template<class T>

PlateBCMap& PlateModel< T >::Impl::getBCMap ( )

inline

Definition at line 555 of file PlateModel_impl.h.

{return _bcMap;}

template<class T>

void PlateModel< T >::Impl::getMoment ( const Mesh &  mesh )

inline

Definition at line 1045 of file PlateModel_impl.h.

References Mesh::getCells(), StorageSite::getCountLevel1(), and sqrt().

  {
    const StorageSite& cells = mesh.getCells();

    const int nCells = cells.getCountLevel1();

    shared_ptr<VectorT3Array> momentPtr(new VectorT3Array(nCells));
    momentPtr->zero();
    _plateFields.moment.addArray(cells,momentPtr);
    VectorT3Array& moment = *momentPtr;

    _deformationGradientModel.compute();

    const VGradArray& wGrad =
      dynamic_cast<const VGradArray&>(_plateFields.deformationGradient[cells]);
    
    const TArray& ym = dynamic_cast<const TArray&>(_plateFields.ym[cells]);
    const TArray& nu = dynamic_cast<const TArray&>(_plateFields.nu[cells]);
    
    const TArray& thickness = dynamic_cast<const TArray&>(_plateFields.thickness[cells]);
    
    TArray& VMStress = dynamic_cast<TArray&>(_plateFields.VMStress[cells]);
    TArray& VMStressOut = dynamic_cast<TArray&>(_plateFields.VMStressOut[cells]);
    VectorT3Array& cellStress =
      dynamic_cast<VectorT3Array&>(_plateFields.stress[cells]);      
    VectorT4Array& devStress =
      dynamic_cast<VectorT4Array&>(_plateFields.devStress[cells]);
    VectorT4Array& pg =
      dynamic_cast<VectorT4Array&>(_plateFields.plasticStrain[cells]);    

    VectorT3Array& plasticMoment =
      dynamic_cast<VectorT3Array&>(_plateFields.plasticMoment[cells]);

    VectorT3Array& cellStrain =
      dynamic_cast<VectorT3Array&>(_plateFields.strain[cells]);
    
    const T onethird(1.0/3.0);
    const T one(1.0);
    const T two(2.0);
    const T three(3.0);
    const T six(6.0);
    const T twelve(12.0);
    const T zero(0.0);

    for(int n=0; n<nCells; n++)
    {
        T cellD = ym[n]*pow(thickness[n],three)/(twelve*(one - pow(nu[n],two)));
        const VGradType& wg = wGrad[n];
        VectorT3 stress;
        
        moment[n][0] = cellD*(wg[0][0]+nu[n]*wg[1][1]);
        moment[n][1] = cellD*(wg[1][1]+nu[n]*wg[0][0]);
        moment[n][2] = cellD*((one-nu[n])/two)*(wg[1][0]+wg[0][1]);
        if(_options.creep)
          moment[n] = moment[n] - plasticMoment[n];
        
        cellStress[n][0] = six*moment[n][0]/pow(thickness[n],two);
        cellStress[n][1] = six*moment[n][1]/pow(thickness[n],two);
        cellStress[n][2] = six*moment[n][2]/pow(thickness[n],two);
       
        cellStrain[n][0] = (thickness[n]/two)*wg[0][0];
        cellStrain[n][1] = (thickness[n]/two)*wg[1][1];
        cellStrain[n][2] = (thickness[n]/two)*(wg[1][0]+wg[0][1]);

        T cellE = ym[n]/(one - pow(nu[n],two));

        int nn = n*(_options.nz+1);
        for(int k=0; k<=_options.nz; k++)
        {
            T zz(thickness[n]*(T(k)-T(_options.nz)/T(2))/T(_options.nz));
            stress[0] = (twelve*zz/pow(thickness[n],three))*cellD*(wg[0][0]+nu[n]*wg[1][1])-
              cellE*(pg[nn+k][0]+nu[n]*pg[nn+k][1]);
            stress[1] = (twelve*zz/pow(thickness[n],three))*cellD*(wg[1][1]+nu[n]*wg[0][0])-
              cellE*(pg[nn+k][1]+nu[n]*pg[nn+k][0]);
            stress[2] = (twelve*zz/pow(thickness[n],three))*cellD*((one-nu[n])/two)*(wg[1][0]+wg[0][1])-
              cellE*(one-nu[n])*(pg[nn+k][3]);
            
            T trace = stress[0] + stress[1];
            devStress[nn+k][0] = stress[0];
            devStress[nn+k][1] = stress[1];
            devStress[nn+k][2] = zero;
            devStress[nn+k][3] = stress[2];
            devStress[nn+k][0] = devStress[nn+k][0] - onethird*trace;
            devStress[nn+k][1] = devStress[nn+k][1] - onethird*trace;
            devStress[nn+k][2] = devStress[nn+k][2] - onethird*trace;
            VMStress[nn+k] = sqrt(pow(stress[0],2.0) + pow(stress[1],2.0) -
                                  stress[0]*stress[1] + three*pow(stress[2],2.0));
            //VMStress[nn+k] = sqrt(pow(stress[0],2.0) + three*pow(stress[2],2.0));
        }
        VMStressOut[n] = VMStress[nn+_options.nz];
      }
  }

template<class T>

PlateModelOptions<T>& PlateModel< T >::Impl::getOptions ( )

inline

Definition at line 557 of file PlateModel_impl.h.

{return _options;}

template<class T>

map<string,shared_ptr<ArrayBase> >& PlateModel< T >::Impl::getPersistenceData ( )

inline

Definition at line 1139 of file PlateModel_impl.h.

References Model::_persistenceData, and Array< T >::zero().

  {
    _persistenceData.clear();
    
    Array<int>* niterArray = new Array<int>(1);
    (*niterArray)[0] = _niters;
    _persistenceData["niters"]=shared_ptr<ArrayBase>(niterArray);

    if (_initialDeformationNorm)
    {
      _persistenceData["initialDeformationNorm"] =
        _initialDeformationNorm->getArrayPtr(_plateFields.deformation);
    }
    else
    {
      Array<Vector<T,3> >* xArray = new Array<Vector<T,3> >(1);
      xArray->zero();
      _persistenceData["initialDeformationNorm"]=shared_ptr<ArrayBase>(xArray);
    }
    return _persistenceData;
  }

template<class T>

PlateVCMap& PlateModel< T >::Impl::getVCMap ( )

inline

Definition at line 556 of file PlateModel_impl.h.

{return _vcMap;}

template<class T>

void PlateModel< T >::Impl::init ( )

inline

Definition at line 376 of file PlateModel_impl.h.

References Model::_meshes, Mesh::getAllFaceGroups(), Mesh::getCells(), StorageSite::getCount(), StorageSite::getCountLevel1(), Mesh::getID(), FaceGroup::groupType, and FaceGroup::site.

  {
    const int numMeshes = _meshes.size();
    for (int n=0; n<numMeshes; n++)
    {
        const Mesh& mesh = *_meshes[n];

        const PlateVC<T>& vc = *_vcMap[mesh.getID()];
            
        const StorageSite& cells = mesh.getCells();
        //        const StorageSite& faces = mesh.getFaces();

        shared_ptr<VectorT3Array> sCell(new VectorT3Array(cells.getCountLevel1()));

        VectorT3 initialDeformation;
        initialDeformation[0] = _options["initialXRotation"];
        initialDeformation[1] = _options["initialYRotation"];
        initialDeformation[2] = _options["initialZDeformation"];
        *sCell = initialDeformation;

        _plateFields.deformation.addArray(cells,sCell);
        _plateFields.deformation.syncLocal();

        if (_options.transient)
        {
            _plateFields.volume0.addArray(cells,
                                         dynamic_pointer_cast<ArrayBase>
                                         (_geomFields.volume[cells].newCopy()));
            _plateFields.volume0.syncLocal();   
                                         
            _plateFields.deformationN1.addArray(cells,
                                            dynamic_pointer_cast<ArrayBase>(sCell->newCopy()));
            _plateFields.deformationN1.syncLocal();
                                                    
            _plateFields.deformationN2.addArray(cells,
                                            dynamic_pointer_cast<ArrayBase>(sCell->newCopy()));
            _plateFields.deformationN2.syncLocal();
                                                                
            if (_options.timeDiscretizationOrder > 1)
              _plateFields.deformationN3.addArray(cells,
                                            dynamic_pointer_cast<ArrayBase>(sCell->newCopy()));
              _plateFields.deformationN3.syncLocal();                                       
                                            
            if(_options.variableTimeStep)
            {
                _options.timeStepN1 = _options["timeStep"];
                _options.timeStepN2 = _options["timeStep"];
            }
        }
        
        shared_ptr<VectorT3Array> stressField(new VectorT3Array(cells.getCountLevel1()));
        stressField->zero();
        _plateFields.stress.addArray(cells,stressField);
        _plateFields.stress.syncLocal();
                
        shared_ptr<VectorT4Array> devStressField(new VectorT4Array((cells.getCountLevel1())*(_options.nz+1)));
        devStressField->zero();
        _plateFields.devStress.addArray(cells,devStressField);
        _plateFields.devStress.syncLocal();

        shared_ptr<TArray> VMStressField(new TArray((cells.getCountLevel1())*(_options.nz+1)));
        VMStressField->zero();
        _plateFields.VMStress.addArray(cells,VMStressField);
        _plateFields.VMStress.syncLocal();

        shared_ptr<TArray> VMStressOutField(new TArray(cells.getCountLevel1()));
        VMStressOutField->zero();
        _plateFields.VMStressOut.addArray(cells,VMStressOutField);
        _plateFields.VMStressOut.syncLocal();

        shared_ptr<VectorT3Array> strainField(new VectorT3Array(cells.getCountLevel1()));
        strainField->zero();
        _plateFields.strain.addArray(cells,strainField);
        _plateFields.strain.syncLocal();        

        shared_ptr<VectorT4Array> plasticStrainField(new VectorT4Array((cells.getCountLevel1())*(_options.nz+1)));
        plasticStrainField->zero();
        _plateFields.plasticStrain.addArray(cells,plasticStrainField);
        _plateFields.plasticStrain.syncLocal();

        shared_ptr<VectorT3Array> plasticStrainOutField(new VectorT3Array(cells.getCountLevel1()));
        plasticStrainOutField->zero();
        _plateFields.plasticStrainOut.addArray(cells,plasticStrainOutField);
        _plateFields.plasticStrainOut.syncLocal();

        shared_ptr<VectorT4Array> plasticStrainN1Field(new VectorT4Array((cells.getCountLevel1())*(_options.nz+1)));
        plasticStrainN1Field->zero();
        _plateFields.plasticStrainN1.addArray(cells,plasticStrainN1Field);
        _plateFields.plasticStrainN1.syncLocal();

        shared_ptr<VectorT3Array> plasticMomentField(new VectorT3Array(cells.getCountLevel1()));
        plasticMomentField->zero();
        _plateFields.plasticMoment.addArray(cells,plasticMomentField); 
        _plateFields.plasticMoment.syncLocal();

        shared_ptr<TArray> rhoCell(new TArray(cells.getCountLevel1()));
        *rhoCell = vc["density"];
        _plateFields.density.addArray(cells,rhoCell);
        _plateFields.density.syncLocal();

        shared_ptr<TArray> ymCell(new TArray(cells.getCountLevel1()));
        *ymCell = vc["ym"];
        _plateFields.ym.addArray(cells,ymCell);
        _plateFields.ym.syncLocal();

        shared_ptr<TArray> nuCell(new TArray(cells.getCountLevel1()));
        *nuCell = vc["nu"];
        _plateFields.nu.addArray(cells,nuCell);
        _plateFields.nu.syncLocal();

        shared_ptr<TArray> forceCell(new TArray(cells.getCountLevel1()));
        forceCell->zero();
        _plateFields.force.addArray(cells,forceCell);
        _plateFields.force.syncLocal();

        shared_ptr<TArray> thicknessCell(new TArray(cells.getCountLevel1()));
        thicknessCell->zero();
        _plateFields.thickness.addArray(cells,thicknessCell);
        _plateFields.thickness.syncLocal();

        shared_ptr<TArray> accelerationCell(new TArray(cells.getCountLevel1()));
        accelerationCell->zero();
        _plateFields.acceleration.addArray(cells,accelerationCell);
        _plateFields.acceleration.syncLocal();
        
        shared_ptr<VectorT3Array> velCell(new VectorT3Array(cells.getCountLevel1()));
        velCell->zero();
        _plateFields.velocity.addArray(cells,velCell);
        _plateFields.velocity.syncLocal();

        //initial temparature gradient array
        shared_ptr<VGradArray> rCell(new VGradArray(cells.getCountLevel1()));
        VGradType residualStress;
        residualStress[0][0] = _options["residualStressXX"];
        residualStress[0][1] = _options["residualStressXY"];
        residualStress[0][2] = _options["residualStressXZ"];
        residualStress[1][0] = _options["residualStressXY"];
        residualStress[1][1] = _options["residualStressYY"];
        residualStress[1][2] = _options["residualStressYZ"];
        residualStress[2][0] = _options["residualStressXZ"];
        residualStress[2][1] = _options["residualStressYZ"];
        residualStress[2][2] = _options["residualStressZZ"];
        *rCell = residualStress;
        _plateFields.residualStress.addArray(cells,rCell);
        _plateFields.residualStress.syncLocal();

        // compute values of deformation flux

        // store deformation flux at interfaces
        /*
        foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
        {
            const FaceGroup& fg = *fgPtr;
            const StorageSite& faces = fg.site;
            shared_ptr<VectorT3Array> deformationFlux(new VectorT3Array(faces.getCount()));
            deformationFlux->zero();
            _plateFields.deformationFlux.addArray(faces,deformationFlux);
        }
        */

        // store deformation flux at boundary faces
        foreach(const FaceGroupPtr fgPtr, mesh.getAllFaceGroups())
        {
            const FaceGroup& fg = *fgPtr;
            const StorageSite& faces = fg.site;
            shared_ptr<VectorT3Array> deformationFlux(new VectorT3Array(faces.getCount()));
            deformationFlux->zero();
            if (fg.groupType != "interior")
            {
                _plateFields.deformationFlux.addArray(faces,deformationFlux);
            }
        }

    }

    _niters  =0;
    _initialDeformationNorm = MFRPtr();
  }

template<class T>

void PlateModel< T >::Impl::initDeformationLinearization ( LinearSystem &  ls )

inline

Definition at line 598 of file PlateModel_impl.h.

References Model::_meshes, MultiField::addArray(), MultiFieldMatrix::addMatrix(), Mesh::getCellCells2(), Mesh::getCells(), LinearSystem::getMatrix(), and LinearSystem::getX().

  {
    const int numMeshes = _meshes.size();
    for (int n=0; n<numMeshes; n++)
    {
        const Mesh& mesh = *_meshes[n];

        const StorageSite& cells = mesh.getCells();
        MultiField::ArrayIndex wIndex(&_plateFields.deformation,&cells);

        ls.getX().addArray(wIndex,_plateFields.deformation.getArrayPtr(cells));

        const CRConnectivity& cellCells2 = mesh.getCellCells2();
        
        shared_ptr<Matrix> m(new CRMatrix<DiagTensorT3,DiagTensorT3,VectorT3>(cellCells2));

        ls.getMatrix().addMatrix(wIndex,wIndex,m);

    }
  }

template<class T>

void PlateModel< T >::Impl::linearizeDeformation ( LinearSystem &  ls )

inline

Definition at line 619 of file PlateModel_impl.h.

References Model::_meshes, PlateBCS< X, Diag, OffDiag >::applyCantileverBC(), PlateBCS< X, Diag, OffDiag >::applyDirichletBC(), PlateBCS< X, Diag, OffDiag >::applyNeumannBC(), PlateBC< T >::bcType, Mesh::getAllFaceGroups(), LinearSystem::getB(), Mesh::getCells(), StorageSite::getCount(), LinearSystem::getMatrix(), MultiFieldMatrix::getMatrix(), PlateModel< T >::getMoment(), FloatVarDict< T >::getVal(), LinearSystem::getX(), FaceGroup::groupType, FaceGroup::id, Linearizer::linearize(), CRMatrix< T_Diag, T_OffDiag, X >::setDirichlet(), and FaceGroup::site.

  {
    _deformationGradientModel.compute();
    DiscrList discretizations;

    const Mesh& mesh = *_meshes[0];
    if(!_options.constForce)
      getMoment(mesh);
    shared_ptr<Discretization>
      sd(new PlateSourceDiscretization<T,DiagTensorT3,DiagTensorT3>
         (_meshes,_geomFields,
          _plateFields.deformation,
          _plateFields.ym,
          _plateFields.nu,
          _plateFields.deformationGradient,
          _plateFields.residualStress,
          _plateFields.thickness,
          _plateFields.force,
          _options.scf,
          _plateFields.devStress,
          _plateFields.VMStress,
          _plateFields.plasticStrain,
          _plateFields.plasticStrainOut,
          _plateFields.plasticStrainN1,
          _plateFields.plasticMoment,
          _options.A,
          _options.B,
          _options.m,
          _options.n,
          _options.Sy0,
          _options.nz,
          _options["timeStep"],
          _options.creepModel,
          _options.creep));
      
    //    shared_ptr<Discretization>
    //  bfd(new SourceDiscretization<VectorT3>
    //     (_meshes,_geomFields,
    //    _structureFields.deformation,
    //    _structureFields.bodyForce));

    discretizations.push_back(sd);
    //discretizations.push_back(bfd);
        
    if (_options.transient)
    {
        shared_ptr<Discretization>
          td(new TimeDerivativePlateDiscretization
             <VectorT3,DiagTensorT3,DiagTensorT3>
             (_meshes,_geomFields,
              _plateFields.deformation,
              _plateFields.deformationN1,
              _plateFields.deformationN2,
              _plateFields.deformationN3,
              _plateFields.density,
              _plateFields.thickness,
              _plateFields.acceleration,
              _plateFields.volume0,
              _options.variableTimeStep,
              _options["timeStep"],
              _options.timeStepN1,
              _options.timeStepN2));
        
        discretizations.push_back(td);
    }
   
    /*
    shared_ptr<Discretization>
      ibm(new GenericIBDiscretization<VectorT3,DiagTensorT3,T>
          (_meshes,_geomFields,_plateFields.deformation));
      
    discretizations.push_back(ibm);
    */
    Linearizer linearizer;

    linearizer.linearize(discretizations,_meshes,ls.getMatrix(),
                         ls.getX(), ls.getB());
    bool allNeumann = true;
    const int numMeshes = _meshes.size();
    for (int n=0; n<numMeshes; n++)
    {
        const Mesh& mesh = *_meshes[n];

        foreach(const FaceGroupPtr fgPtr, mesh.getAllFaceGroups())
        {
            const FaceGroup& fg = *fgPtr;
            if (fg.groupType != "interior")
            {
            const StorageSite& faces = fg.site;
            const int nFaces = faces.getCount();
            const TArray& faceAreaMag = 
              dynamic_cast<const TArray&>(_geomFields.areaMag[faces]);
            const PlateBC<T>& bc = *_bcMap[fg.id];
            

            PlateBCS<VectorT3,DiagTensorT3,DiagTensorT3> gbc(faces,mesh,
                                                    _geomFields,
                                                    _plateFields.deformation,
                                                    ls.getMatrix(), ls.getX(), ls.getB());
            
            VectorT3 fluxB(NumTypeTraits<VectorT3>::getZero());
            VectorT3 zeroFlux(NumTypeTraits<VectorT3>::getZero());
            if (bc.bcType == "Clamped")
            {
                FloatValEvaluator<VectorT3>
                  bDeformation(bc.getVal("specifiedXRotation"),
                               bc.getVal("specifiedYRotation"),
                               bc.getVal("specifiedZDeformation"),
                               faces);
                fluxB = gbc.applyDirichletBC(bDeformation);
                
                allNeumann = false;
            }
            else if (bc.bcType == "Symmetry")
            {
                gbc.applySymmetryBC();
                allNeumann = false;
            }
            else if (bc.bcType == "ZeroDerivative")
            {
                gbc.applyZeroDerivativeBC();
            }
            else if (bc.bcType == "SpecifiedTraction")
            {
                for(int f=0; f<nFaces; f++)
                {
                    
                    fluxB += gbc.applyNeumannBC(f,zeroFlux);
        
                }
            }
            else if (bc.bcType == "SpecifiedForce")
            {
                FloatValEvaluator<VectorT3>
                  bForce(bc.getVal("specifiedXForce"),
                         bc.getVal("specifiedYForce"),
                         bc.getVal("specifiedZForce"),
                         faces);
                for(int f=0; f<nFaces; f++)
                {

                    fluxB += gbc.applyNeumannBC(f,bForce[f]/faceAreaMag[f]);

                }
            }
            else if (bc.bcType == "SpecifiedShear")
            {
                FloatValEvaluator<VectorT3>
                  bShear(bc.getVal("specifiedXShear"),
                         bc.getVal("specifiedYShear"),
                         bc.getVal("specifiedZShear"),
                         faces);
                for(int f=0; f<nFaces; f++)
                {
                        
                    fluxB += gbc.applyCantileverBC(f,bShear[f]);

                }
            }
            else if (bc.bcType != "Interface")
              throw CException(bc.bcType + " not implemented for StructureModel");
            //cout << "force sum for " << fg.id  << " = " << fluxB << endl;  
            }
        }
        /*
        foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
        {
            const FaceGroup& fg = *fgPtr;
            const StorageSite& faces = fg.site;
            GenericBCS<VectorT3,DiagTensorT3,DiagTensorT3> gbc(faces,mesh,
                                                    _geomFields,
                                                    _structureFields.deformation,
                                                    _structureFields.deformationFlux,
                                                    ls.getMatrix(), ls.getX(), ls.getB());

            gbc.applyInterfaceBC();
        }
        */
    }
        
        
#ifdef FVM_PARALLEL
     int count = 1;
     int allNeumannInt = int( allNeumann);
     MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &allNeumannInt, count, MPI::INT, MPI::PROD);
     allNeumann = bool(allNeumannInt);
#endif
        
        
                
    if(allNeumann && !_options.transient)
    {
        const Mesh& mesh = *_meshes[0];
        const StorageSite& cells = mesh.getCells();
        MultiField& b = ls.getB();
        MultiFieldMatrix& matrix = ls.getMatrix();
        MultiField::ArrayIndex wIndex(&_plateFields.deformation,&cells);
        VectorT3Array& rCell = dynamic_cast<VectorT3Array&>(b[wIndex]);
        VectorT3Array& w = dynamic_cast<VectorT3Array&>
          (_plateFields.deformation[cells]);
        VVMatrix& vvMatrix =
          dynamic_cast<VVMatrix&>(matrix.getMatrix(wIndex,wIndex));
        rCell[0] = T(0);
        w[0] = T(0);
        vvMatrix.setDirichlet(0);
    }

#if 0
    shared_ptr<Discretization>
      ud(new Underrelaxer<VectorT3,DiagTensorT3,DiagTensorT3>
         (_meshes,_structureFields.deformation,
          _options["deformationURF"]));
    
    DiscrList discretizations2;
    discretizations2.push_back(ud);

    linearizer.linearize(discretizations2,_meshes,ls.getMatrix(),
                         ls.getX(), ls.getB());
#endif

  }

template<class T>

void PlateModel< T >::Impl::postPlateSolve ( LinearSystem &  ls )

inline

Definition at line 980 of file PlateModel_impl.h.

References Model::_meshes, Mesh::getCells(), StorageSite::getCountLevel1(), and LinearSystem::getDelta().

  {
    MultiField& sField = ls.getDelta();

    const int numMeshes = _meshes.size();
    for(int n=0;n<numMeshes;n++)
    {
        const Mesh& mesh = *_meshes[n];

        const StorageSite& cells = mesh.getCells();

        MultiField::ArrayIndex sIndex(&_plateFields.deformation,&cells);
        VectorT3Array& w = dynamic_cast<VectorT3Array&>
          (_plateFields.deformation[cells]);
        const VectorT3Array& ww = dynamic_cast<const VectorT3Array&>
          (sField[sIndex]);
        //const T deformationURF(_options["deformationURF"]);

        const int nCells = cells.getCountLevel1();
        for(int c=0;c<nCells;c++)
        {
            w[c] += ww[c];
        }
    }
  }

template<class T>

void PlateModel< T >::Impl::printBCs ( )

inline

Definition at line 1032 of file PlateModel_impl.h.

  {
    foreach(typename PlateBCMap::value_type& pos, _bcMap)
    {
        cout << "Face Group " << pos.first << ":" << endl;
        cout << "    bc type " << pos.second->bcType << endl;
        foreach(typename PlateBC<T>::value_type& vp, *pos.second)
        {
            cout << "   " << vp.first << " "  << vp.second.constant <<  endl;
        }
    }
  }

template<class T>

void PlateModel< T >::Impl::restart ( )

inline

Definition at line 1161 of file PlateModel_impl.h.

References Model::_persistenceData.

  {
    if (_persistenceData.find("niters") != _persistenceData.end())
    {
        shared_ptr<ArrayBase> rp = _persistenceData["niters"];
        ArrayBase& r = *rp;
        Array<int>& niterArray = dynamic_cast<Array<int>& >(r);
        _niters = niterArray[0];
    }
    if (_persistenceData.find("initialDeformationNorm") != _persistenceData.end())
    {
        shared_ptr<ArrayBase>  r = _persistenceData["initialDeformationNorm"];
        _initialDeformationNorm = MFRPtr(new MultiFieldReduction());
        _initialDeformationNorm->addArray(_plateFields.deformation,r);
    }
  }

template<class T>

MFRPtr PlateModel< T >::Impl::solveDeformation ( )

inline

Definition at line 842 of file PlateModel_impl.h.

References LinearSystem::initAssembly(), LinearSystem::initSolve(), LinearSystem::postSolve(), and LinearSystem::updateSolution().

  {
    LinearSystem ls;
    initDeformationLinearization(ls);
    ls.initAssembly();
    linearizeDeformation(ls);
    ls.initSolve();    
    
    //AMG solver(ls);
    MFRPtr rNorm = _options.getDeformationLinearSolver().solve(ls);
    if (!_initialDeformationNorm) _initialDeformationNorm = rNorm;
        
    _options.getDeformationLinearSolver().cleanup();
    ls.postSolve();
    ls.updateSolution();
    //postPlateSolve(ls);
    if (_options.transient){
      calculatePlateVelocity();
      calculatePlateAcceleration();
    }
    return rNorm;
  }

template<class T>

void PlateModel< T >::Impl::updateTime ( )

inline

Definition at line 559 of file PlateModel_impl.h.

References Model::_meshes, and Mesh::getCells().

  {
    const int numMeshes = _meshes.size();
    for (int n=0; n<numMeshes; n++)
    {
        const Mesh& mesh = *_meshes[n];

        const StorageSite& cells = mesh.getCells();
        VectorT3Array& w =
          dynamic_cast<VectorT3Array&>(_plateFields.deformation[cells]);
        VectorT3Array& wN1 =
          dynamic_cast<VectorT3Array&>(_plateFields.deformationN1[cells]);
        VectorT3Array& wN2 =
          dynamic_cast<VectorT3Array&>(_plateFields.deformationN2[cells]);
        if (_options.timeDiscretizationOrder > 1)
        {
            VectorT3Array& wN3 =
              dynamic_cast<VectorT3Array&>(_plateFields.deformationN3[cells]);
            wN3 = wN2;
        }
        wN2 = wN1;
        wN1 = w;
        if(_options.variableTimeStep)
        {
            if (_options.timeDiscretizationOrder > 1)
            {
                _options.timeStepN2 = _options.timeStepN1;
            }
            _options.timeStepN1 = _options["timeStep"];
        }
        VectorT4Array& pS =
          dynamic_cast<VectorT4Array&>(_plateFields.plasticStrain[cells]);
        VectorT4Array& pSN1 =
          dynamic_cast<VectorT4Array&>(_plateFields.plasticStrainN1[cells]);
        pSN1 = pS;
    }
  }

Member Data Documentation

template<class T>

PlateBCMap PlateModel< T >::Impl::_bcMap

private

Definition at line 1281 of file PlateModel_impl.h.

template<class T>

GradientModel<VectorT3> PlateModel< T >::Impl::_deformationGradientModel

private

Definition at line 1285 of file PlateModel_impl.h.

template<class T>

const GeomFields& PlateModel< T >::Impl::_geomFields

private

Definition at line 1278 of file PlateModel_impl.h.

template<class T>

MFRPtr PlateModel< T >::Impl::_initialDeformationNorm

private

Definition at line 1287 of file PlateModel_impl.h.

template<class T>

const MeshList PlateModel< T >::Impl::_meshes

private

Definition at line 1277 of file PlateModel_impl.h.

template<class T>

int PlateModel< T >::Impl::_niters

private

Definition at line 1288 of file PlateModel_impl.h.

template<class T>

PlateModelOptions<T> PlateModel< T >::Impl::_options

private

Definition at line 1284 of file PlateModel_impl.h.

template<class T>

map<string,shared_ptr<ArrayBase> > PlateModel< T >::Impl::_persistenceData

private

Definition at line 1290 of file PlateModel_impl.h.

template<class T>

PlateFields& PlateModel< T >::Impl::_plateFields

private

Definition at line 1279 of file PlateModel_impl.h.

template<class T>

PlateVCMap PlateModel< T >::Impl::_vcMap

private

Definition at line 1282 of file PlateModel_impl.h.

---

The documentation for this class was generated from the following file:

• PlateModel_impl.h