composite_properties_lop.hh

#ifndef COMP_PROP_LOP_HH
#define COMP_PROP_LOP_HH
 
#include <dune/geometry/type.hh>
#include <dune/geometry/referenceelements.hh>
#include <dune/geometry/quadraturerules.hh>
 
namespace Dune {
  namespace PDELab {
 
    template<typename GV, typename MODEL>
    class composite_properties :
        public NumericalJacobianApplyVolume<composite_properties<GV,MODEL> >,
        public FullVolumePattern,
        public LocalOperatorDefaultFlags,
        public InstationaryLocalOperatorDefaultMethods<double>,
        public NumericalJacobianVolume<composite_properties<GV,MODEL> >
    {
    public:
      // pattern assembly flags
      enum { doPatternVolume = false };
 
      // residual assembly flags
      enum { doAlphaVolume  = true };
 
      //Constructor
      composite_properties (const GV& gv_, const MODEL& model_) : model(model_), gv(gv_) {}
      //
      // alpha_volume for getStress
      //
      // Volume Integral Depending on Test and Ansatz Functions
      template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
        void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
        {
          //Unwrap stress
          const auto& lfsv1 = lfsv.template child<0>();
          const auto& lfsv2 = lfsv.template child<1>();
 
          int id = gv.indexSet().index(eg.entity());
 
          int materialType = model.getMaterialTypeFromElement(id);
 
          int orientation = model.getOrientation(id);
 
          r.accumulate(lfsv1, 0, materialType);
          r.accumulate(lfsv2, 0, orientation);
        }
    private:
      const MODEL& model;
      const GV& gv;
    };
 
  }
}
 
#endif