plot_properties.hh

#ifndef PLOT_STRESSES_HH
#define PLOT_STRESSES_HH
 
#include "composite_properties_lop.hh"
 
namespace Dune{
    namespace Composites{
 
        template<class MODEL, typename V, typename GV1, typename GFS, typename MBE>
            void plotProperties(MODEL& model, const GV1& gv1, const GFS& gfs,  MBE& mbe){
                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);
 
                typedef typename GV::Grid::ctype e_ctype;
                typedef Dune::PDELab::QkDGLocalFiniteElementMap<e_ctype,double,0,3> FEM_properties;
                FEM_properties fem_properties;
 
                typedef Dune::PDELab::GridFunctionSpace<GV, FEM_properties, Dune::PDELab::NoConstraints, Scalar_VectorBackend> SCALAR_P1GFS;
                SCALAR_P1GFS materialType(gv,fem_properties); materialType.name("materialType");
                SCALAR_P1GFS orientation(gv,fem_properties); orientation.name("orientation");
 
                typedef Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::fixed> VectorBackend;
 
                typedef Dune::PDELab::PowerGridFunctionSpace <SCALAR_P1GFS, 2, VectorBackend,Dune::PDELab::EntityBlockedOrderingTag> P1GFS;
 
                P1GFS p1gfs(materialType,orientation);
 
                typedef Dune::PDELab::EmptyTransformation NoTrafo;
 
                Dune::PDELab::composite_properties<GV,MODEL> lopProperties(gv,model);
 
                typedef Dune::PDELab::GridOperator<P1GFS,P1GFS,Dune::PDELab::composite_properties<GV,MODEL>,MBE,RF,RF,RF,NoTrafo,NoTrafo> GO;
                GO go(p1gfs,p1gfs,lopProperties,mbe);
 
                typedef typename GO::Traits::Range V1;
                V1 properties(p1gfs,0.0), dummy(p1gfs,0.0); 
                go.residual(dummy,properties);
 
                Dune::SubsamplingVTKWriter<GV1> vtkwriter(gv1,Dune::refinementLevels(0));
                Dune::PDELab::addSolutionToVTKWriter(vtkwriter,p1gfs,properties);
                vtkwriter.write(model.vtk_properties_output,Dune::VTK::appendedraw);
            }
 
 
                template<class MODEL, typename V, typename GV1, typename GFS, typename MBE>
            void plotPropertiesUG(MODEL& model, const GV1& gv1, const GFS& gfs,  MBE& mbe){
                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);
 
                typedef typename GV::Grid::ctype e_ctype;
                typedef Dune::PDELab::PkLocalFiniteElementMap<GV,e_ctype,double,1> FEM_properties;
                FEM_properties fem_properties(gv);
 
                typedef Dune::PDELab::GridFunctionSpace<GV, FEM_properties, Dune::PDELab::NoConstraints, Scalar_VectorBackend> SCALAR_P1GFS;
                SCALAR_P1GFS materialType(gv,fem_properties); materialType.name("materialType");
                SCALAR_P1GFS orientation(gv,fem_properties); orientation.name("orientation");
 
                typedef Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::fixed> VectorBackend;
 
                typedef Dune::PDELab::PowerGridFunctionSpace <SCALAR_P1GFS, 2, VectorBackend,Dune::PDELab::EntityBlockedOrderingTag> P1GFS;
 
                // Define operator to count elements containing a degree of freeedom
                Dune::PDELab::countElem lopCnt;
                SCALAR_P1GFS cntpgfs(gv,fem_properties);
                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);
 
                V x0(gfs,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 plot properties
                P1GFS p1gfs(materialType,orientation);
 
                typedef Dune::PDELab::EmptyTransformation NoTrafo;
 
                Dune::PDELab::composite_propertiesUG<GV,MODEL, DGF> lopProperties(gv,model,dgfCnt);
 
                typedef Dune::PDELab::GridOperator<P1GFS,P1GFS,Dune::PDELab::composite_propertiesUG<GV,MODEL,DGF>,MBE,RF,RF,RF,NoTrafo,NoTrafo> GO;
                GO go(p1gfs,p1gfs,lopProperties,mbe);
 
                typedef typename GO::Traits::Range V1;
                V1 properties(p1gfs,0.0), dummy(p1gfs,0.0); 
                go.residual(dummy,properties);
 
                Dune::SubsamplingVTKWriter<GV1> vtkwriter(gv1,Dune::refinementLevels(0));
                Dune::PDELab::addSolutionToVTKWriter(vtkwriter,p1gfs,properties);
                vtkwriter.write(model.vtk_properties_output,Dune::VTK::appendedraw);
            }
 
    }
}
 
#endif