composite_properties_lop.hh

// -*- tab-width: 2; indent-tabs-mode: nil -*-
// vi: set et ts=2 sw=2 sts=2:
#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, typename DGF>
    class composite_propertiesUG :
      public NumericalJacobianApplyVolume<composite_propertiesUG<GV,MODEL,DGF> >,
      public FullVolumePattern,
      public LocalOperatorDefaultFlags,
      public InstationaryLocalOperatorDefaultMethods<double>,
      public NumericalJacobianVolume<composite_propertiesUG<GV,MODEL,DGF> >
    {
    public:
      // pattern assembly flags
      enum { doPatternVolume = false };
 
      // residual assembly flags
      enum { doAlphaVolume  = true };
 
      //Constructor
      composite_propertiesUG (const GV& gv_, const MODEL& model_, const DGF& dgf_) : model(model_), gv(gv_), dgf(dgf_) {}
      //
      // 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
      {
        [[maybe_unused]] typedef typename LFSV::Traits::SizeType size_type;
 
        //Unwrap stress
        const auto& lfsv1 = lfsv.template child<0>();
        const auto& lfsv2 = lfsv.template child<1>();
 
        int id = gv.indexSet().index(eg.entity());
 
        // size_type nodes_per_element = lfsv1.size();
        // size_type dof = lfsv1.size();
        int materialType, orientation;
        materialType = model.getMaterialTypeFromElement(id);
        orientation = model.getOrientation(id);
        // select quadrature rule
        assert(eg.geometry().corners() == static_cast<int>(lfsv1.size()));
        assert(eg.geometry().corners() == static_cast<int>(lfsv2.size()));
        for (int iter = 0; iter < eg.geometry().corners(); iter++)
        {
          Dune::FieldVector<double,3> pt = eg.geometry().local(eg.geometry().corner(iter));
          Dune::FieldVector<double,1> num_elem(1.0);
          dgf.evaluate(eg.entity(),pt,num_elem);
          r.accumulate(lfsv1, iter, materialType/num_elem[0]);
          r.accumulate(lfsv2, iter, orientation/num_elem[0]);
        }
 
      }
    private:
      const MODEL& model;
      const GV& gv;
      const DGF& dgf;
    };
 
    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, orientation;
        materialType = model.getMaterialTypeFromElement(id);
        orientation = model.getOrientation(id);
        r.accumulate(lfsv1, 0, materialType);
        r.accumulate(lfsv2, 0, orientation);
      }
    private:
      const MODEL& model;
      const GV& gv;
    };
  }
}
 
#endif