computeStresses.hh

#ifndef COMPUTE_STRESSES_HH
#define COMPUTE_STRESSES_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<Dune::PDELab::istl::Blocking::none,1> 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,6> VectorBackend;
    typedef Dune::PDELab::CompositeGridFunctionSpace <VectorBackend,Dune::PDELab::EntityBlockedOrderingTag, 
        SCALAR_P1GFS, SCALAR_P1GFS, SCALAR_P1GFS, SCALAR_P1GFS, SCALAR_P1GFS, SCALAR_P1GFS> 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,0);
                Dune::PDELab::addSolutionToVTKWriter(vtkwriter,p1gfs,stress);
                vtkwriter.write(model.vtk_stress_output,Dune::VTK::appendedraw);
    }
}
 
#endif