d3f_debug_linker.h

Go to the documentation of this file.

/* SPDX-FileCopyrightText: 2025 Gesellschaft fuer Anlagen- und Reaktorsicherheit gGmbH
 * SPDX-License-Identifier: EUPL-1.2
 * SPDX-FileContributor: Andreas Vogel, Arne Naegel
 * SPDX-FileContributor: Goethe Universität Frankfurt
 * SPDX-FileType: SOURCE
 * 
 * This file is part of d3f++.
 * d3f++ is an extension for UG4. Licensing information and citation requirements of UG4 are provided in LICENSES/UG4-LGPL_2.1
 */
 
#ifndef __H__UG__LIB_DISC__SPATIAL_DISC__D3F_LINKER__
#define __H__UG__LIB_DISC__SPATIAL_DISC__D3F_LINKER__
 
 
 
#include "lib_disc/common/geometry_util.h"
#include "lib_disc/spatial_disc/user_data/std_user_data.h"
#include "lib_disc/spatial_disc/user_data/linker/linker.h"
#include "common/math/math_vector_matrix/math_vector.h"
 
 
 
 
namespace ug{
 
 
// Peclet number linker
 
template <int dim>
class DFPecletNumber
  : public StdDataLinker< DFPecletNumber<dim>, number, dim>
{
    typedef StdDataLinker< DFPecletNumber<dim>, number, dim> base_type;
 
  public:
    DFPecletNumber() :
      m_spDispersion(NULL), /*m_spDDispersion(NULL),*/
      m_spVelocity(NULL), /*m_spDVelocity(NULL),*/
      m_dLengthScale(1.0)
    {
    //  this linker has exactly two inputs
      this->set_num_input(2);
    }
 
 
    inline void evaluate (number& value,
                          const MathVector<dim>& globIP,
                          number time, int si) const
    {
      MathVector<dim> q, Dq;
      MathMatrix<dim,dim> D;
 
      (*m_spVelocity)(q, globIP, time, si);
      (*m_spDispersion)(D, globIP, time, si);
 
      number L = m_dLengthScale;
      number qnorm = VecTwoNorm(q);
 
      // Compute Diffusion Tensor times Gradient
      MatVecMult(Dq, D, q);
      const number qtDq = VecDot(Dq, q);
 
      value = (qnorm*qnorm*qnorm*L)/qtDq;
 
    //  UG_ASSERT(value >= -1e-8, "Error");
    //  std::cout << VecTwoNorm(velocity) << " " << MinAbsEigenvalue(dispersion) <<  " "<< value <<std::endl;
      UG_LOG("evaluate1:" << value << std::endl);
    }
 
    template <int refDim>
    inline void evaluate(number vValue[],
                         const MathVector<dim> vGlobIP[],
                         number time, int si,
                         GridObject* elem,
                         const MathVector<dim> vCornerCoords[],
                         const MathVector<refDim> vLocIP[],
                         const size_t nip,
                         LocalVector* u,
                         const MathMatrix<refDim, dim>* vJT = NULL) const
    {
      ReferenceObjectID roid=elem->reference_object_id();
      number size = ElementSize<dim>(roid, vCornerCoords);
      number hL = pow(size, 1.0/dim)*m_dLengthScale;
 
      std::vector<MathVector<dim> > vVelocity(nip);
      std::vector<MathMatrix<dim,dim> > vDispersion(nip);
 
      (*m_spVelocity)(&vVelocity[0], vGlobIP, time, si,
                elem, vCornerCoords, vLocIP, nip, u, vJT);
      (*m_spDispersion)(&vDispersion[0], vGlobIP, time, si,
              elem, vCornerCoords, vLocIP, nip, u, vJT);
 
      for(size_t ip = 0; ip < nip; ++ip)
      {
 
        const number qnorm_ip = VecTwoNorm(vVelocity[ip]);
 
        MathVector<dim> Dvel_ip;
        MatVecMult(Dvel_ip, vDispersion[ip], vVelocity[ip]);
        const number qtDq_ip = VecDot(Dvel_ip, vVelocity[ip]);
 
        vValue[ip] = qnorm_ip*qnorm_ip*qnorm_ip*hL;
        vValue[ip] /= qtDq_ip;
 
        //UG_ASSERT(vValue[ip] >= -1e-8, "Error");
        //std::cout << "ip=" << ip <<" : "<<VecTwoNorm(vVelocity[ip]) << " " << MinAbsEigenvalue(vDispersion[ip]) << "="<< vValue[ip] << std::endl;
        // UG_LOG("evaluate2:" << vValue[ip] << std::endl);
      }
    }
 
    template <int refDim>
    void eval_and_deriv(number vValue[],
                        const MathVector<dim> vGlobIP[],
                        number time, int si,
                        GridObject* elem,
                        const MathVector<dim> vCornerCoords[],
                        const MathVector<refDim> vLocIP[],
                        const size_t nip,
                        LocalVector* u,
                        bool bDeriv,
                        int s,
                        std::vector<std::vector<number > > vvvDeriv[]//,
                        //const MathVector<refDim>* vJT = NULL
                        ) const
    {
 
      //  Compute the derivatives at all ips     //
 
    //  check if something to do
      if(!bDeriv || this->zero_derivative()) return;
 
      UG_ASSERT(0, "eval_and_deriv not implemented.");
    //  clear all derivative values
      this->set_zero(vvvDeriv, nip);
 
    }
 
  public:
    void set_dispersion(SmartPtr<CplUserData<MathMatrix<dim,dim>, dim> > data)
    {
      m_spDispersion = data;
      //m_spDDispersion = data.template cast_dynamic<DependentUserData<MathMatrix<dim,dim>, dim> >();
      base_type::set_input(_DISPERSION_, data, data);
    }
 
    void set_dispersion(number val)
    { set_dispersion(make_sp(new ConstUserMatrix<dim>(val))); }
 
    void set_velocity(SmartPtr<CplUserData<MathVector<dim>, dim> > data)
    {
      m_spVelocity = data;
      // m_spDVelocity = data.template cast_dynamic<DependentUserData<MathVector<dim>, dim> >();
      base_type::set_input(_VELOCITY_, data, data);
    }
 
    void set_length_scale(number l)
    { m_dLengthScale=l; }
 
  protected:
    enum {_DISPERSION_ = 0,  _VELOCITY_ = 1};
    SmartPtr<CplUserData<MathMatrix<dim,dim>, dim> > m_spDispersion;
    // SmartPtr<DependentUserData<MathMatrix<dim,dim>, dim> > m_spDDispersion;
 
    SmartPtr<CplUserData<MathVector<dim>, dim> > m_spVelocity;
    // SmartPtr<DependentUserData<MathVector<dim>, dim> > m_spDVelocity;
 
    number m_dLengthScale;
 
};
 
 
template <int dim>
class DFAnisotropy
    : public StdDataLinker< DFAnisotropy<dim>, number, dim>
{
    typedef StdDataLinker< DFAnisotropy<dim>, number, dim> base_type;
 
    public:
      DFAnisotropy() : m_spDispersion(NULL)/*, m_spDDispersion(NULL)*/
      {
          this->set_num_input(1);
      }
 
 
      inline void evaluate (number& value,
                            const MathVector<dim>& globIP,
                            number time, int si) const
      {
        MathMatrix<dim,dim> dispersion;
        (*m_spDispersion)(dispersion, globIP, time, si);
 
        value = MinAbsEigenvalue(dispersion);
        value /= MaxAbsEigenvalue(dispersion);
        UG_LOG("aniso_evaluate1:" << value << std::endl);
      }
 
      template <int refDim>
      inline void evaluate(number vValue[],
                           const MathVector<dim> vGlobIP[],
                           number time, int si,
                           GridObject* elem,
                           const MathVector<dim> vCornerCoords[],
                           const MathVector<refDim> vLocIP[],
                           const size_t nip,
                           LocalVector* u,
                           const MathMatrix<refDim, dim>* vJT = NULL) const
      {
        std::vector<MathMatrix<dim,dim> > vDispersion(nip);
        (*m_spDispersion)(&vDispersion[0], vGlobIP, time, si,
                elem, vCornerCoords, vLocIP, nip, u, vJT);
 
        for(size_t ip = 0; ip < nip; ++ip)
        {
          vValue[ip]  = MinAbsEigenvalue(vDispersion[ip]);
          vValue[ip] /= MaxAbsEigenvalue(vDispersion[ip]);
          UG_LOG("aniso_evaluate2:" << vValue[ip] << std::endl);
        }
 
      }
 
      template <int refDim>
      void eval_and_deriv(number vValue[],
                          const MathVector<dim> vGlobIP[],
                          number time, int si,
                          GridObject* elem,
                          const MathVector<dim> vCornerCoords[],
                          const MathVector<refDim> vLocIP[],
                          const size_t nip,
                          LocalVector* u,
                          bool bDeriv,
                          int s,
                          std::vector<std::vector<number > > vvvDeriv[]//,
                          //const MathVector<refDim>* vJT = NULL
                          ) const
      {
 
        //  Compute the derivatives at all ips     //
 
      //  check if something to do
        if(!bDeriv || this->zero_derivative()) return;
 
        UG_ASSERT(0, "eval_and_deriv not implemented.");
      //  clear all derivative values
        this->set_zero(vvvDeriv, nip);
 
      }
 
    public:
      void set_dispersion(SmartPtr<CplUserData<MathMatrix<dim,dim>, dim> > data)
      {
        m_spDispersion = data;
        //m_spDDispersion = data.template cast_dynamic<DependentUserData<MathMatrix<dim,dim>, dim> >();
        base_type::set_input(_DISPERSION_, data, data);
      }
 
      void set_dispersion(number val)
      { set_dispersion(make_sp(new ConstUserMatrix<dim>(val))); }
 
    protected:
      enum {_DISPERSION_ = 0};
      SmartPtr<CplUserData<MathMatrix<dim,dim>, dim> > m_spDispersion;
      //SmartPtr<DependentUserData<MathMatrix<dim,dim>, dim> > m_spDDispersion;
 
  };
 
template <int dim>
class ElementExtensions
  //: public StdDataLinker< ElementExtensions<dim>, MathVector<dim>, dim>
  : public StdDependentUserData<ElementExtensions<dim>, MathVector<dim>, dim>
{
  public:
    typedef StdDependentUserData<ElementExtensions<dim>, MathVector<dim>, dim> base_type;
    typedef MathVector<dim> data_type;
 
  public:
    ElementExtensions()
    {}
/*
    virtual void operator() (MathVector<dim>& value,
                 const MathVector<dim>& globIP,
                 number time, int si) const
    { UG_THROW("global operator() not implemented."); }
 
 
    virtual void operator()(MathVector<dim> vValue[],
                const MathVector<dim> vGlobIP[],
                number time, int si, const size_t nip) const
    { UG_THROW("global operator() not implemented."); }
 
*/
    /*inline void evaluate (number& value,
                              const MathVector<dim>& globIP,
                              number time, int si) const
    {}*/
/*
    template <int refDim>
    inline void evaluate(MathVector<dim> vValue[],
                         const MathVector<dim> vGlobIP[],
                         number time, int si,
                         GridObject* elem,
                         const MathVector<dim> vCornerCoords[],
                         const MathVector<refDim> vLocIP[],
                         const size_t nip,
                         LocalVector* u,
                         const MathMatrix<refDim, dim>* vJT = NULL) const
    {
      ReferenceObjectID roid=elem->reference_object_id();
      MathVector<dim> ext;
      ElementExtensionsSq<dim>(roid, ext, vCornerCoords);
      for(size_t ip = 0; ip < nip; ++ip)
        vValue[ip] = ext;
 
    }
    */
    template <int refDim>
    void eval_and_deriv(data_type vValue[],
                              const MathVector<dim> vGlobIP[],
                              number time, int si,
                              GridObject* elem,
                              const MathVector<dim> vCornerCoords[],
                              const MathVector<refDim> vLocIP[],
                              const size_t nip,
                              LocalVector* u,
                              bool bDeriv,
                              int s,
                              std::vector<std::vector<data_type > > vvvDeriv[],
                              const MathMatrix<refDim, dim>* vJT = NULL
                              ) const
    {
      ReferenceObjectID roid=elem->reference_object_id();
      MathVector<dim> ext;
      ElementExtensionsSq<dim>(roid, ext, vCornerCoords);
      for(size_t ip = 0; ip < nip; ++ip)
      {
              vValue[ip] = ext;
      }
 
 
      //  Compute the derivatives at all ips     //
 
 
      //  check if something to do
      if(!bDeriv || this->zero_derivative()) return;
 
      UG_ASSERT(0, "eval_and_deriv not implemented.");
          //  clear all derivative values
            //this->set_zero(vvvDeriv, nip);
 
    }
 
/*
    virtual void compute(LocalVector* u, GridObject* elem,
                     const MathVector<dim> vCornerCoords[], bool bDeriv = false)
    {
      ReferenceObjectID roid=elem->reference_object_id();
      MathVector<dim> ext;
      ElementExtensionsSq<dim>(roid, ext, vCornerCoords);
 
      for(size_t s = 0; s < this->num_series(); ++s)
        for(size_t ip = 0; ip < this->num_ip(s); ++ip)
          this->value(s,ip) = ext;
      UG_ASSERT(0, "TODO: Implementation as 'UserData' is missing.");
    }
 
    virtual void compute(LocalVectorTimeSeries* u, GridObject* elem,
                 const MathVector<dim> vCornerCoords[], bool bDeriv = false)
    {
      ReferenceObjectID roid=elem->reference_object_id();
      MathVector<dim> ext;
      ElementExtensionsSq<dim>(roid, ext, vCornerCoords);
 
      UG_ASSERT(0, "TODO: Implementation as 'UserData' is missing.");
 
    }*/
 
    virtual bool constant() const {return false;}
 
    virtual bool requires_grid_fct() const {return true;}
 
    virtual bool continuous() const {return false;}
 
 
};
 
} // end namespace ug
 
#endif /* __H__UG__LIB_DISC__SPATIAL_DISC__DARCY_VELOCITY_LINKER__ */