DGtal  1.5.beta
DirichletConditions.h
1
17 #pragma once
18
31 #if defined(DirichletConditions_RECURSES)
32 #error Recursive header files inclusion detected in DirichletConditions.h
33 #else // defined(DirichletConditions_RECURSES)
35 #define DirichletConditions_RECURSES
36
37 #if !defined DirichletConditions_h
39 #define DirichletConditions_h
40
42 // Inclusions
43 #include <iostream>
44 #include <vector>
45 #include "DGtal/base/Common.h"
46 #include <DGtal/math/linalg/CDynamicMatrix.h>
47 #include <DGtal/math/linalg/CDynamicVector.h>
48 #include <DGtal/math/linalg/CLinearAlgebra.h>
50
51 namespace DGtal
52 {
53
55  // template class DirichletConditions
96  template < typename TLinearAlgebraBackend >
98  {
99  public:
100  typedef TLinearAlgebraBackend LinearAlgebraBackend;
101
103  typedef typename LinearAlgebraBackend::DenseVector::Scalar Scalar;
104  typedef typename LinearAlgebraBackend::DenseVector DenseVector;
105  typedef typename LinearAlgebraBackend::IntegerVector IntegerVector;
108  typedef typename LinearAlgebraBackend::Triplet Triplet;
109
115
129  static
131  const IntegerVector& p )
132  {
133  ASSERT( A.cols() == A.rows() );
134  ASSERT( p.rows() == A.rows() );
135  const auto n = p.rows();
136  std::vector< Index > relabeling( n );
137  Index j = 0;
138  for ( Index i = 0; i < p.rows(); i++ )
139  relabeling[ i ] = ( p[ i ] == 0 ) ? j++ : n;
140  // Building matrix
141  SparseMatrix Ap( j, j );
142  std::vector< Triplet > triplets;
143  for ( int k = 0; k < A.outerSize(); ++k )
144  for ( typename SparseMatrix::InnerIterator it( A, k ); it; ++it )
145  {
146  if ( ( relabeling[ it.row() ] != n ) && ( relabeling[ it.col() ] != n ) )
147  triplets.push_back( { relabeling[ it.row() ], relabeling[ it.col() ],
148  it.value() } );
149  }
150  Ap.setFromTriplets( triplets.cbegin(), triplets.cend() );
151  return Ap;
152  }
153
169  static
171  const DenseVector& b,
172  const IntegerVector& p,
173  const DenseVector& u )
174  {
175  ASSERT( A.cols() == A.rows() );
176  ASSERT( p.rows() == A.rows() );
177  const auto n = p.rows();
178  DenseVector up = p.array().template cast<double>() * u.array();
179  DenseVector tmp = b.array() - (A * up).array();
180  std::vector< Index > relabeling( n );
181  Index j = 0;
182  for ( Index i = 0; i < p.rows(); i++ )
183  relabeling[ i ] = ( p[ i ] == 0 ) ? j++ : n;
184  DenseVector bp( j );
185  for ( Index i = 0; i < p.rows(); i++ )
186  if ( p[ i ] == 0 ) bp[ relabeling[ i ] ] = tmp[ i ];
187  return bp;
188  }
189
202  static
204  const IntegerVector& p,
205  const DenseVector& u )
206  {
207  DenseVector x( p.rows() );
208  Index j = 0;
209  for ( Index i = 0; i < p.rows(); i++ )
210  x[ i ] = ( p[ i ] == 0 ) ? xd[ j++ ] : u[ i ];
211  return x;
212  }
213
222  static
224  {
225  return IntegerVector::Zero( b.rows() );
226  }
227  };
228
229 } // namespace DGtal
230
231
233 // Includes inline functions.
234
235 // //
237
238 #endif // !defined DirichletConditions_h
239
240 #undef DirichletConditions_RECURSES
241 #endif // else defined(DirichletConditions_RECURSES)
242
Aim: A helper class to solve a system with Dirichlet boundary conditions.
BOOST_CONCEPT_ASSERT((concepts::CDynamicMatrix< SparseMatrix >))
LinearAlgebraBackend::DenseMatrix DenseMatrix
LinearAlgebraBackend::Triplet Triplet
LinearAlgebraBackend::IntegerVector IntegerVector
BOOST_CONCEPT_ASSERT((concepts::CDynamicVector< DenseVector >))
static DenseVector dirichletVector(const SparseMatrix &A, const DenseVector &b, const IntegerVector &p, const DenseVector &u)
static SparseMatrix dirichletOperator(const SparseMatrix &A, const IntegerVector &p)
BOOST_CONCEPT_ASSERT((concepts::CDynamicMatrix< DenseMatrix >))
BOOST_CONCEPT_ASSERT((concepts::CLinearAlgebra< DenseVector, SparseMatrix >))
LinearAlgebraBackend::DenseVector::Index Index
TLinearAlgebraBackend LinearAlgebraBackend
static IntegerVector nullBoundaryVector(const DenseVector &b)
LinearAlgebraBackend::DenseVector::Scalar Scalar
static DenseVector dirichletSolution(const DenseVector &xd, const IntegerVector &p, const DenseVector &u)
BOOST_CONCEPT_ASSERT((concepts::CLinearAlgebra< DenseVector, DenseMatrix >))
LinearAlgebraBackend::SparseMatrix SparseMatrix
LinearAlgebraBackend::DenseVector DenseVector
SMesh::Index Index
DGtal is the top-level namespace which contains all DGtal functions and types.
Aim: Represent any dynamic sized matrix having sparse or dense representation.
Aim: Represent any dynamic sized column vector having sparse or dense representation.
Aim: Check right multiplication between matrix and vector and internal matrix multiplication....
EigenLinearAlgebraBackend::SparseMatrix SparseMatrix
EigenLinearAlgebraBackend::DenseMatrix DenseMatrix