symmetricmatrix.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright © DUNE-FUFEM Project contributors, see file AUTHORS.md
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
 
#ifndef DUNE_FUFEM_SYMMETRICMATRIX_HH
#define DUNE_FUFEM_SYMMETRICMATRIX_HH
 
#include <dune/common/fvector.hh>
 
namespace Dune {
namespace Fufem {
 
template <class T, int N>
class SymmetricMatrix
{
public:
 
  using field_type = T;
 
  using Data = Dune::FieldVector<T,N*(N+1)/2>;
 
  SymmetricMatrix()
  {}
 
 
  template <class Iterator>
  SymmetricMatrix(Iterator it)
  {
    for (size_t i=0; i<data_.size(); ++i, ++it)
      data_[i] = *it;
  }
 
  static constexpr SymmetricMatrix<T,N> identityMatrix()
  {
    SymmetricMatrix<T,N> id;
    for( size_t i=0; i<N; ++i)
      id.setEntry(i,i,1);
 
    return id;
  }
 
  static constexpr Data traceMap()
  {
    // trace(A) = <Id,A>, therefore return data of Id
    // use isometry!
    return identityMatrix().data();
  }
 
  SymmetricMatrix<T,N>& operator=(const T& s)
  {
    data_ = s;
    return *this;
  }
 
  SymmetricMatrix<T,N>& operator*=(const T& s)
  {
    data_ *= s;
    return *this;
  }
 
  void setEntry(int i, int j, const T& entry)
  {
    if( i > j )
      data_[i*(i+1)/2+j] = std::sqrt(2.) * entry;
    else if ( i < j )
      data_[j*(j+1)/2+i] = std::sqrt(2.) * entry;
    else // ( i == j )
      data_[i*(i+1)/2+i] = entry;
  }
 
  T operator() (int i, int j) const
  {
    if( i > j )
      return 1./std::sqrt(2.) * data_[i*(i+1)/2+j];
    else if ( i < j )
      return 1./std::sqrt(2.) * data_[j*(j+1)/2+i];
    else // ( i == j )
      return data_[i*(i+1)/2+i];
  }
 
  T energyScalarProduct (const FieldVector<T,N>& v1, const FieldVector<T,N>& v2) const
  {
    T result = 0;
    for (size_t i=0; i<N; i++)
      for (size_t j=0; j<=i; j++)
          result += (1-0.5*(i==j)) * operator()(i,j) * (v1[i] * v2[j] + v1[j] * v2[i]);
    return result;
  }
 
  void axpy(const T& a, const SymmetricMatrix<T,N>& other)
  {
    data_.axpy(a,other.data_);
  }
 
  T frobenius_norm2() const
  {
    return data_.two_norm2();
  }
 
  T frobenius_norm() const
  {
    return std::sqrt(frobenius_norm2());
  }
 
  const Data& data() const
  {
    return data_;
  }
 
  Data& data()
  {
    return data_;
  }
 
  static constexpr size_t dataSize()
  {
    return Data::size();
  }
 
private:
  Data data_;
};
 
} // namespace Fufem
} // namespace Dune
 
#endif