computeStresses.hh

#ifndef COMPUTE_STRESSES_HH
#define COMPUTE_STRESSES_HH
 
#include<dune/pdelab/finiteelementmap/qkdg.hh>
#include "getStress.hh"
 
template<class MODEL, typename V, typename GV1, typename GFS, typename MBE>
void calculateStresses(MODEL& model, V& x0, const GV1& gv1, const GFS& gfs,  MBE& mbe, int dofel = 60){
    using Dune::PDELab::Backend::native;
 
    typedef Dune::PDELab::ISTL::VectorBackend<> Scalar_VectorBackend;
    typedef double RF;
    std::string vtk_output = config.get<std::string>("ModelName","TestModel");
    Dune::Timer timer;
 
    using GV = Dune::PDELab::AllEntitySet<GV1>; 
    auto gv = GV(gv1);
 
    // Initially we consider only a DG stress plot
   //Function spaces for Stresses
    typedef typename GV::Grid::ctype e_ctype;
    typedef Dune::PDELab::QkDGLocalFiniteElementMap<e_ctype,double,0,3> FEM_stress;
        FEM_stress fem_stress;
 
    typedef Dune::PDELab::GridFunctionSpace<GV, FEM_stress, Dune::PDELab::NoConstraints, Scalar_VectorBackend> SCALAR_P1GFS;
    SCALAR_P1GFS p1gfs_00(gv,fem_stress); p1gfs_00.name("stress00");  SCALAR_P1GFS p1gfs_11(gv,fem_stress); p1gfs_11.name("stress11");
    SCALAR_P1GFS p1gfs_22(gv,fem_stress); p1gfs_22.name("stress22");  SCALAR_P1GFS p1gfs_12(gv,fem_stress); p1gfs_12.name("stress12");
    SCALAR_P1GFS p1gfs_02(gv,fem_stress); p1gfs_02.name("stress02");  SCALAR_P1GFS p1gfs_01(gv,fem_stress); p1gfs_01.name("stress01");
 
    typedef Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::fixed> VectorBackend;
    typedef Dune::PDELab::PowerGridFunctionSpace <SCALAR_P1GFS,6,VectorBackend,Dune::PDELab::EntityBlockedOrderingTag> P1GFS;
 
    P1GFS p1gfs(p1gfs_00,p1gfs_11,p1gfs_22,p1gfs_12,p1gfs_02,p1gfs_01);
 
    //Set up local operator to calculate stresses
    Dune::PDELab::getStress<GV,MODEL,60> lopStress(gv,model);
 
    typedef Dune::PDELab::EmptyTransformation NoTrafo;
    typedef Dune::PDELab::GridOperator<GFS,P1GFS,Dune::PDELab::getStress<GV,MODEL,60>,MBE,RF,RF,RF,NoTrafo,NoTrafo> GOStress;
        GOStress goStress(gfs,p1gfs,lopStress,mbe);
 
    typedef typename GOStress::Traits::Range V1;
        V1 stress(p1gfs,0.0); 
        goStress.residual(x0,stress);
 
    // Record Stresses to MODEL class to be used later if required
    model.sizeStressVector(native(stress).size()); // Make sure stress contain in model class is the right size
    for (unsigned int i = 0; i < native(stress).size(); i++){
        model.setStress(i,native(stress)[i]);
    }
 
    // Save solution in vtk format for paraview
    if(model.stress_Plot){
                Dune::SubsamplingVTKWriter<GV1> vtkwriter(gv1,Dune::refinementLevels(0));
                Dune::PDELab::addSolutionToVTKWriter(vtkwriter,p1gfs,stress);
                vtkwriter.write(model.vtk_stress_output,model.stress_Plot_ascii ? Dune::VTK::ascii : Dune::VTK::appendedraw);
    }
}
 
 
template<class MODEL, typename V, typename GV1, typename GFS, typename MBE>
void calculateStressesUG(MODEL& model, V& x0, const GV1& gv1, const GFS& gfs,  MBE& mbe, int dofel = 12){
    using Dune::PDELab::Backend::native;
    using GV = Dune::PDELab::AllEntitySet<GV1>; 
    auto gv = GV(gv1);
    typedef double RF;
 
    typedef typename GV::Grid::ctype e_ctype;
    typedef Dune::PDELab::PkLocalFiniteElementMap<GV,e_ctype,double,1> FEM_stress;
    // Construct FEM Space
    FEM_stress fem_stress(gv);
    typedef Dune::PDELab::ISTL::VectorBackend<> Scalar_VectorBackend;
 
    typedef Dune::PDELab::GridFunctionSpace<GV, FEM_stress, Dune::PDELab::NoConstraints, Scalar_VectorBackend> SCALAR_P1GFS;
    SCALAR_P1GFS p1gfs_00(gv,fem_stress); p1gfs_00.name("stress00");  SCALAR_P1GFS p1gfs_11(gv,fem_stress); p1gfs_11.name("stress11");
    SCALAR_P1GFS p1gfs_22(gv,fem_stress); p1gfs_22.name("stress22");  SCALAR_P1GFS p1gfs_12(gv,fem_stress); p1gfs_12.name("stress12");
    SCALAR_P1GFS p1gfs_02(gv,fem_stress); p1gfs_02.name("stress02");  SCALAR_P1GFS p1gfs_01(gv,fem_stress); p1gfs_01.name("stress01");
    typedef Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::fixed> VectorBackend;
    typedef Dune::PDELab::PowerGridFunctionSpace <SCALAR_P1GFS,6,VectorBackend,Dune::PDELab::EntityBlockedOrderingTag> P1GFS;
 
    P1GFS p1gfs(p1gfs_00,p1gfs_11,p1gfs_22,p1gfs_12,p1gfs_02,p1gfs_01);
 
    Dune::PDELab::countElem lopCnt;
    SCALAR_P1GFS cntpgfs(gv,fem_stress);
    typedef Dune::PDELab::EmptyTransformation NoTrafo;
    typedef Dune::PDELab::GridOperator<GFS,SCALAR_P1GFS,Dune::PDELab::countElem,MBE,RF,RF,RF,NoTrafo,NoTrafo> GOCnt;
    GOCnt goCnt(gfs,cntpgfs,lopCnt,mbe);
 
    typedef typename GOCnt::Traits::Range Vcnt;
    Vcnt ElemCnt(cntpgfs,0.0);
 
    goCnt.residual(x0,ElemCnt);
 
    //Turn counter into discrete grid function
    typedef Dune::PDELab::DiscreteGridFunction<SCALAR_P1GFS,Vcnt> DGF;
    DGF dgfCnt(cntpgfs, ElemCnt);
 
    //Set up local operator to calculate stresses
    double maxDisp = 0.0;
    Dune::PDELab::getStressUG<GV,MODEL,DGF,12> lopStress(gv,model,dgfCnt,maxDisp);
    typedef Dune::PDELab::GridOperator<GFS,P1GFS,Dune::PDELab::getStressUG<GV,MODEL,DGF,12>,MBE,RF,RF,RF,NoTrafo,NoTrafo> GOStress;
    GOStress goStress(gfs,p1gfs,lopStress,mbe);
 
    typedef typename GOStress::Traits::Range V1;
    V1 stress(p1gfs,0.0);
    goStress.residual(x0,stress);
 
    // Record Stresses to MODEL class to be used later if required
    model.sizeStressVector(native(stress).size()); // Make sure stress contain in model class is the right size
    for (unsigned int i = 0; i < native(stress).size(); i++){
        model.setStress(i,native(stress)[i]);
    }
 
    // Save solution in vtk format for paraview
    if(model.stress_Plot){
                Dune::SubsamplingVTKWriter<GV1> vtkwriter(gv1,Dune::refinementLevels(0));
                Dune::PDELab::addSolutionToVTKWriter(vtkwriter,p1gfs,stress);
                vtkwriter.write(model.vtk_stress_output,model.stress_Plot_ascii ? Dune::VTK::ascii : Dune::VTK::appendedraw);
    }
}
 
#endif