df_transform.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: 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__DF_TRANSFORM_H_
#define __H__UG__DF_TRANSFORM_H_
 
#include <vector>
 
 
 
namespace ug{
namespace df{
 
 
 
 
 
 
template<class TAlgebra>
class DFTransform
{
 
public:
 
 
  typedef typename TAlgebra::matrix_type matrix_type;
  typedef typename TAlgebra::vector_type vector_type;
 
protected:
 
  enum DFComponents {CONCENTRATION=0, PRESSURE=1};
 
 
  static int transform_operator(const std::vector<double> &sigma, const matrix_type &A, matrix_type &Anew)
  {
 
      // std::cout << "DFLeftTransformMatrix....";
 
      // copy
      Anew.set_as_copy_of(A);
 
      // eliminate off-diag entries of diagonal element
      const size_t N = A.num_rows();
      for(size_t i=0; i < N; i++)
      {
        // eliminate off-diag entries of diagonal element
        // "eliminate" i-th row of A
        const typename TAlgebra::matrix_type::row_iterator rowEnd = Anew.end_row(i);
        for(typename TAlgebra::matrix_type::row_iterator ij_iter = Anew.begin_row(i);
            ij_iter != rowEnd; ++ij_iter)
        {
          typename TAlgebra::matrix_type::value_type &Aij = ij_iter.value();
          //std::cout << Aij << "->" << sigma[i];
 
 
          for (size_t beta=0; beta<GetCols(Aij); ++beta)
          {
            BlockRef(Aij, CONCENTRATION, beta) -= sigma[i] * BlockRef(Aij, PRESSURE, beta);
          }
 
          //std::cout << "->" << Aij << std::endl;
        }
      }
 
      return 0;
  };
 
 
  static int transform_rhs(const std::vector<double> &sigma, const vector_type &rhs, vector_type &rhs_new)
  {
      // std::cout << "DFLeftTransformRHS....";
 
      // copy
      rhs_new=rhs;
 
      // transform
      const size_t N = rhs_new.size();
      for(size_t i=1; i < N; i++)
      {
        // eliminate off-diag entries of diagonal element
        BlockRef(rhs_new[i], CONCENTRATION) -= sigma[i] * BlockRef(rhs_new[i], PRESSURE);
        //std::cout << BlockRef(rhs[i], CONCENTRATION) <<  "->"<< BlockRef(rhs_new[i], CONCENTRATION) << std::endl;
      }
      return 0;
  };
 
};
 
template<class TAlgebra>
class DFAlgebraicLeftTransform : public DFTransform<TAlgebra>
{
 
 
public:
  typedef typename TAlgebra::matrix_type matrix_type;
  typedef typename TAlgebra::vector_type vector_type;
  DFAlgebraicLeftTransform() : sigma(){}
 
  // called once
  void init(const matrix_type &A, const vector_type *sol)
  {
    const size_t N = A.num_rows();
    sigma.resize(N);
 
    for(size_t i=1; i < N; i++)
    {
      const typename TAlgebra::matrix_type::value_type &Aii=A(i,i);
 
      sigma[i] = BlockRef(Aii, DFTransform<TAlgebra>::CONCENTRATION, DFTransform<TAlgebra>::PRESSURE) /
                BlockRef(Aii, DFTransform<TAlgebra>::PRESSURE, DFTransform<TAlgebra>::PRESSURE);
    }
  }
 
  int transform_operator(const matrix_type &A, matrix_type &Anew)
  { return DFTransform<TAlgebra>::transform_operator(sigma, A, Anew); }
 
  int transform_rhs(const vector_type &rhs, vector_type &rhs_new)
  { return DFTransform<TAlgebra>::transform_rhs(sigma, rhs, rhs_new); }
 
protected:
  std::vector<double> sigma;
};
 
 
template<class TAlgebra>
class DFLeftTransform : public DFTransform<TAlgebra>
{
 
 
public:
  typedef typename TAlgebra::matrix_type matrix_type;
  typedef typename TAlgebra::vector_type vector_type;
 
  DFLeftTransform() : sigma(){}
 
  void init(const matrix_type &A, const vector_type *sol)
  {
    const size_t N = A.num_rows();
    sigma.resize(N);
 
    for(size_t i=1; i < N; i++)
    {
      const typename TAlgebra::matrix_type::value_type &Aii=A(i,i);
 
      sigma[i] = BlockRef(Aii, DFTransform<TAlgebra>::CONCENTRATION, DFTransform<TAlgebra>::PRESSURE) /
                BlockRef(Aii, DFTransform<TAlgebra>::PRESSURE, DFTransform<TAlgebra>::PRESSURE);
    }
  }
 
  int transform_operator(const matrix_type &A, matrix_type &Anew)
  { return DFTransform<TAlgebra>::transform_operator(sigma, A, Anew); }
 
  int transform_rhs(const vector_type &rhs, vector_type &rhs_new)
  { return DFTransform<TAlgebra>::transform_rhs(sigma, rhs, rhs_new); }
 
protected:
  std::vector<double> sigma;
};
 
 
 
} // namespace df
} // namespace ug
#endif /* DF_PRECOND_H_ */