DirichletConditions.h

 
#pragma once
 
#if defined(DirichletConditions_RECURSES)
#error Recursive header files inclusion detected in DirichletConditions.h
#else // defined(DirichletConditions_RECURSES)
#define DirichletConditions_RECURSES
 
#if !defined DirichletConditions_h
#define DirichletConditions_h
 
// Inclusions
#include <iostream>
#include <vector>
#include "DGtal/base/Common.h"
#include <DGtal/math/linalg/CDynamicMatrix.h>
#include <DGtal/math/linalg/CDynamicVector.h>
#include <DGtal/math/linalg/CLinearAlgebra.h>
 
namespace DGtal
{
 
  // template class DirichletConditions
  template < typename TLinearAlgebraBackend >
  class DirichletConditions
  {
  public:
    typedef TLinearAlgebraBackend LinearAlgebraBackend;
 
    typedef typename LinearAlgebraBackend::DenseVector::Index  Index;
    typedef typename LinearAlgebraBackend::DenseVector::Scalar Scalar;
    typedef typename LinearAlgebraBackend::DenseVector         DenseVector;
    typedef typename LinearAlgebraBackend::IntegerVector       IntegerVector;
    typedef typename LinearAlgebraBackend::DenseMatrix         DenseMatrix;
    typedef typename LinearAlgebraBackend::SparseMatrix        SparseMatrix;
    typedef typename LinearAlgebraBackend::Triplet             Triplet;
 
    BOOST_CONCEPT_ASSERT(( concepts::CDynamicVector<DenseVector> ));
    BOOST_CONCEPT_ASSERT(( concepts::CDynamicMatrix<DenseMatrix> ));
    BOOST_CONCEPT_ASSERT(( concepts::CDynamicMatrix<SparseMatrix> ));
    BOOST_CONCEPT_ASSERT(( concepts::CLinearAlgebra<DenseVector, SparseMatrix> ));
    BOOST_CONCEPT_ASSERT(( concepts::CLinearAlgebra<DenseVector, DenseMatrix> ));
 
    static 
    SparseMatrix dirichletOperator( const SparseMatrix& A,
                                    const IntegerVector& p )
    {
      ASSERT( A.cols() == A.rows() );
      ASSERT( p.rows() == A.rows() );
      const auto n = p.rows();
      std::vector< Index > relabeling( n );
      Index j = 0;
      for ( Index i = 0; i < p.rows(); i++ )
        relabeling[ i ] = ( p[ i ] == 0 ) ? j++ : n;
      // Building matrix
      SparseMatrix Ap( j, j );
      std::vector< Triplet > triplets;
      for ( int k = 0; k < A.outerSize(); ++k )
        for ( typename SparseMatrix::InnerIterator it( A, k ); it; ++it )
          {
            if ( ( relabeling[ it.row() ] != n ) && ( relabeling[ it.col() ] != n ) )
              triplets.push_back( { relabeling[ it.row() ], relabeling[ it.col() ],
                  it.value() } );
          }
      Ap.setFromTriplets( triplets.cbegin(), triplets.cend() );
      return Ap;
    }
 
    static
    DenseVector dirichletVector( const SparseMatrix& A,
                                 const DenseVector& b,
                                 const IntegerVector& p,
                                 const DenseVector& u ) 
    {
      ASSERT( A.cols() == A.rows() );
      ASSERT( p.rows() == A.rows() );
      const auto n = p.rows();
      DenseVector  up = p.array().template cast<double>() * u.array();
      DenseVector tmp = b.array() - (A * up).array();
      std::vector< Index > relabeling( n );
      Index j = 0;
      for ( Index i = 0; i < p.rows(); i++ )
        relabeling[ i ] = ( p[ i ] == 0 ) ? j++ : n;
      DenseVector  bp( j );
      for ( Index i = 0; i < p.rows(); i++ )
        if ( p[ i ] == 0 ) bp[ relabeling[ i ] ] = tmp[ i ];
      return bp;
    }
 
    static
    DenseVector dirichletSolution( const DenseVector& xd,
                                   const IntegerVector& p,
                                   const DenseVector& u ) 
    {
      DenseVector  x( p.rows() );
      Index j = 0;
      for ( Index i = 0; i < p.rows(); i++ )
        x[ i ] = ( p[ i ] == 0 ) ? xd[ j++ ] : u[ i ];
      return x;
    }
 
    static
    IntegerVector nullBoundaryVector( const DenseVector& b )
    {
      return IntegerVector::Zero( b.rows() );
    }
  };
 
} // namespace DGtal
 
 
// Includes inline functions.
 
//                                                                           //
 
#endif // !defined DirichletConditions_h
 
#undef DirichletConditions_RECURSES
#endif // else defined(DirichletConditions_RECURSES)