DGtal  1.3.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)
34 
35 #define DirichletConditions_RECURSES
36 
37 #if !defined DirichletConditions_h
38 
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 
DGtal::DirichletConditions::dirichletVector
static DenseVector dirichletVector(const SparseMatrix &A, const DenseVector &b, const IntegerVector &p, const DenseVector &u)
Definition: DirichletConditions.h:170
DGtal::concepts::CDynamicVector
Aim: Represent any dynamic sized column vector having sparse or dense representation.
Definition: CDynamicVector.h:88
Index
SMesh::Index Index
Definition: fullConvexitySphereGeodesics.cpp:117
DGtal::DirichletConditions::SparseMatrix
LinearAlgebraBackend::SparseMatrix SparseMatrix
Definition: DirichletConditions.h:107
DGtal::DirichletConditions::DenseVector
LinearAlgebraBackend::DenseVector DenseVector
Definition: DirichletConditions.h:104
DGtal::DirichletConditions::dirichletSolution
static DenseVector dirichletSolution(const DenseVector &xd, const IntegerVector &p, const DenseVector &u)
Definition: DirichletConditions.h:203
DGtal::DirichletConditions::dirichletOperator
static SparseMatrix dirichletOperator(const SparseMatrix &A, const IntegerVector &p)
Definition: DirichletConditions.h:130
DGtal::DirichletConditions::Triplet
LinearAlgebraBackend::Triplet Triplet
Definition: DirichletConditions.h:108
SparseMatrix
EigenLinearAlgebraBackend::SparseMatrix SparseMatrix
Definition: testHeatLaplace.cpp:50
DGtal::DirichletConditions::DenseMatrix
LinearAlgebraBackend::DenseMatrix DenseMatrix
Definition: DirichletConditions.h:106
DGtal::DirichletConditions::nullBoundaryVector
static IntegerVector nullBoundaryVector(const DenseVector &b)
Definition: DirichletConditions.h:223
DGtal
DGtal is the top-level namespace which contains all DGtal functions and types.
DGtal::DirichletConditions::BOOST_CONCEPT_ASSERT
BOOST_CONCEPT_ASSERT((concepts::CDynamicVector< DenseVector >))
DGtal::DirichletConditions::IntegerVector
LinearAlgebraBackend::IntegerVector IntegerVector
Definition: DirichletConditions.h:105
DenseMatrix
EigenLinearAlgebraBackend::DenseMatrix DenseMatrix
Definition: testHeatLaplace.cpp:51
DGtal::concepts::CDynamicMatrix
Aim: Represent any dynamic sized matrix having sparse or dense representation.
Definition: CDynamicMatrix.h:88
DGtal::DirichletConditions
Aim: A helper class to solve a system with Dirichlet boundary conditions.
Definition: DirichletConditions.h:97
DGtal::DirichletConditions::Scalar
LinearAlgebraBackend::DenseVector::Scalar Scalar
Definition: DirichletConditions.h:103
DGtal::concepts::CLinearAlgebra
Aim: Check right multiplication between matrix and vector and internal matrix multiplication....
Definition: CLinearAlgebra.h:93
DGtal::DirichletConditions::Index
LinearAlgebraBackend::DenseVector::Index Index
Definition: DirichletConditions.h:102
DGtal::DirichletConditions::LinearAlgebraBackend
TLinearAlgebraBackend LinearAlgebraBackend
Definition: DirichletConditions.h:100