linearStaticDriver.hh

#include <dune/composites/Setup/parallelPartition.hh>
#include <dune/composites/Setup/GridTransformation.hh>
#include <dune/composites/Setup/dirichletBC.hh>
#include <dune/composites/Driver/FEM/Serendipity/serendipityfem.hh>
 
#include "../PostProcessing/computeStresses.hh"
#include "../PostProcessing/plot_properties.hh"
#include "localOperator/linearelasticity.hh"
 
namespace Dune{
    namespace Composites{
 
        template <typename MODEL>
            class linearStaticDriver{
 
                public:
                    typedef Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::none,1> Scalar_VectorBackend;
                    typedef Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::fixed,3> VectorBackend;
                    typedef Dune::PDELab::istl::BCRSMatrixBackend<> MBE;
                    typedef double RF;
 
                    int elem_order;
 
                    Dune::MPIHelper& helper;
 
                    linearStaticDriver(Dune::MPIHelper & helper_) : helper(helper_){ 
                        elem_order = 2; // Default Value which Can be overwritten by setElementOrder
                    };
 
                    void inline setElementOrder(int val = 2){
                        elem_order = val;
                    }
 
                    void inline apply(MODEL& myModel){
                        if (myModel.setUp_Required == true){
                            myModel.LayerCake();
                            myModel.setUpMaterials();
                            myModel.computeElasticTensors();
                        }
 
                        Dune::Timer watch;
                        std::vector<double> times(3);
 
                        // === A: Setup YaspGrid
                        watch.reset();
                        typedef Dune::YaspGrid<3,Dune::TensorProductCoordinates<double,3>> YGRID;
                        YaspPartition<3> yp(helper.size(),myModel.GridPartition());
                        YGRID yaspgrid(myModel.coords,myModel.periodic,myModel.overlap,helper.getCommunicator(),(Dune::YLoadBalance<3>*)&yp);
                        int refinements = myModel.refineBaseGrid;
                        if(refinements > 0) yaspgrid.globalRefine(refinements);
                        typedef YGRID::LeafGridView YGV;
                        const YGV& ygv = yaspgrid.leafGridView();
 
                        if(helper.rank() == 0){
                            std::cout << "Number of elements per processor: " << ygv.size(0) << std::endl;
                            std::cout << "Number of nodes per processor: "    << ygv.size(3) << std::endl;
                        }
 
                        myModel.setPG(ygv); // Loops over elements and assigns a physical group
 
                        // ==================================================================
                        //                         Transform Grid 
                        // ==================================================================
                        typedef Dune::Composites::GridTransformation<3,MODEL> GRID_TRAFO;
                        GRID_TRAFO gTrafo(myModel, helper.rank());
 
                        typedef typename Dune::GeometryGrid<YGRID,GRID_TRAFO> GRID;
                        GRID grid(yaspgrid,gTrafo);
                        if(helper.rank() == 0)
                            std::cout << "Grid transformation complete" << std::endl;
 
                        //Define Grid view
                        typedef typename GRID::LeafGridView GV;
                        const GV& gv = grid.leafGridView();
 
                        if(helper.rank() == 0)
                            std::cout << "Grid view set up" << std::endl;
 
                        times[1] = watch.elapsed();
                        watch.reset();
 
                        typedef Scalar_BC<GV,MODEL,RF> BC;
                        typedef Dune::PDELab::CompositeConstraintsParameters<BC,BC,BC> Constraints;
 
                        // ==================================================================
                        //                         Set up problem
                        // ==================================================================
                        myModel.template computeTensorsOnGrid<GV,YGV>(gv,ygv);
 
                        // Setup initial boundary conditions for each degree of freedom
                        typedef Scalar_BC<GV,MODEL,double> InitialDisp;
                        InitialDisp u1(gv, myModel,0), u2(gv, myModel,1), u3(gv, myModel,2);
 
                        // Wrap scalar boundary conditions in to vector
                        typedef Dune::PDELab::CompositeGridFunction<InitialDisp,InitialDisp,InitialDisp>        InitialSolution;
                        InitialSolution initial_solution(u1,u2,u3);
 
                        // Construct grid function spaces for each degree of freedom
                        typedef Dune::PDELab::OverlappingConformingDirichletConstraints CON;
                        CON con;
 
                        times[2] = watch.elapsed();
                        watch.reset();
 
                        if(elem_order == 2){ // == Element Order
                            const int element_order = 2;
                            const int dofel = 3 * 20;
                            const int non_zeros = 81;
 
                            if(helper.rank() == 0)
                                std::cout << "Piecewise quadratic serendipity elements" << std::endl;
 
                            typedef Dune::PDELab::SerendipityLocalFiniteElementMap<GV,typename GV::Grid::ctype,double,element_order> FEM;
                            FEM fem(gv);
 
                            typedef Dune::PDELab::GridFunctionSpace<GV, FEM, CON, Scalar_VectorBackend> SCALAR_GFS;
                            SCALAR_GFS dispU1(gv,fem,con); dispU1.name("U1");
                            SCALAR_GFS dispU2(gv,fem,con); dispU2.name("U2");
                            SCALAR_GFS dispU3(gv,fem,con); dispU3.name("U3");
 
                            // Note that Vectors are blocked by Dune::PDELab::EntityBlockedOrderingTag
                            typedef Dune::PDELab::CompositeGridFunctionSpace <VectorBackend,Dune::PDELab::EntityBlockedOrderingTag, SCALAR_GFS, SCALAR_GFS, SCALAR_GFS> GFS;
                            const GFS gfs(dispU1, dispU2, dispU3);
 
                            typedef typename GFS::template ConstraintsContainer<RF>::Type C;
 
                            // Make constraints map and initialize it from a function
                            C cg;
                            cg.clear();
 
                            BC U1_cc(gv,myModel,0), U2_cc(gv,myModel,1), U3_cc(gv,myModel,2);
                            Constraints constraints(U1_cc,U2_cc,U3_cc);
 
                            Dune::PDELab::constraints(constraints,gfs,cg);
 
                            MBE mbe(non_zeros); // Maximal number of nonzeroes per row
 
                            // === Construct Linear Operator on FEM Space
                            typedef Dune::PDELab::linearelasticity<GV, MODEL, dofel> LOP;
                            LOP lop(gv, myModel);
 
                            typedef Dune::PDELab::GridOperator<GFS,GFS,LOP,MBE,RF,RF,RF,C,C> GO;
                            GO go(gfs,cg,gfs,cg,lop,mbe);
 
                            // === Make coefficent vector and initialize it from a function
                            typedef Dune::PDELab::Backend::Vector<GFS,double> V;
                            V u(gfs,0.0);
                            Dune::PDELab::interpolate(initial_solution,gfs,u);
 
                            // === Set non constrained dofs to zero
                            Dune::PDELab::set_shifted_dofs(cg,0.0,u);
                            myModel.template solve<GO,V,GFS,C,Constraints,MBE,LOP>(go,u,gfs,cg,constraints,mbe,lop);
 
 
                            Dune::SubsamplingVTKWriter<GV> vtkwriter(gv,0);
                            Dune::PDELab::addSolutionToVTKWriter(vtkwriter,gfs,u);
                            vtkwriter.write(myModel.vtk_displacement_output,Dune::VTK::appendedraw);
 
 
                            calculateStresses<MODEL,V,GV,GFS,MBE>(myModel,u,gv,gfs,mbe);
 
                            plotProperties<MODEL,V,GV,GFS,MBE>(myModel,gv,gfs,mbe);
 
                            myModel.template postprocess<GO,V,GFS,C,MBE,GV>(go,u,gfs,cg,gv,mbe);
                        }
                    } // end apply
 
            };
    }
}