1 /*
2 XLiFE++ is an extended library of finite elements written in C++
3 Copyright (C) 2014 Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin
4
5 This program is free software: you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation, either version 3 of the License, or
8 (at your option) any later version.
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13 You should have received a copy of the GNU General Public License
14 along with this program. If not, see <http://www.gnu.org/licenses/>.
15 */
16
17 /*!
18 \file TermMatrix.hpp
19 \author E. Lunéville
20 \since 23 mar 2012
21 \date 5 oct 2012
22
23 \brief Definition of the xlifepp::TermMatrix class
24
25 xlifepp::TermMatrix class handles the numerical representation of general bilinear forms (xlifepp::BilinearForm)
26 that is the matrix a(vi,wj), i=1,...m; j=1,..,n with (vi)and (wj) the basis functions.
27 It inherits from xlifepp::Term class and it is an end user class.
28
29 A xlifepp::BilinearForm object (see BilinearForm.hpp) is a multiple unknowns/spaces bilinear form
30 that is a list of bilinear forms acting on a couple of (sub)spaces (xlifepp::SuBilinearForm object)
31 indexed by unknown/testfunction . As a consquence, a xlifepp::TermMatrix object is a list of xlifepp::SuTermMatrix objects
32 that are related to xlifepp::SuBilinearForm objects
33
34 \verbatim
35 child
36 Term --------> TermMatrix
37 | BilinearForm
38 | MatrixEntry*
39 | map<uvPair,SuTermMatrix> (uvPair is a pair of unknown, testfunction)
40 | ...
41 \endverbatim
42
43 End user has to instanciate only xlifepp::TermMatrix object, xlifepp::SuTermMatrix is an internal object never managed by end user,
44 even if its bilinear form is really a single unknown/testfunction bilinear form.
45
46 xlifepp::MatrixEntry pointer of xlifepp::TermMatrix is in general null, may be not in case of multiple uvPair bilinear form.
47 xlifepp::TermMatrix being seen as a block matrix, its matrix values are given by xlifepp::MatrixEntry of each xlifepp::SuTermMatrix.
48 In some circonstances (direct solvers), xlifepp::MatrixEntry pointer of xlifepp::TermMatrix may be allocated to store xlifepp::TermMatrix
49 as a contiguous matrix without block structure.
50 Be care, both representations may exist at the same time but with two times memory consuming !!!
51
52 Note : if required it is possible to go to scalar representation of matrix by using scalar_entries_p
53 and renumbering vectors cdofs_u (column) and cdofs_v (row).
54 This representation has no interest if all the unknowns involved in xlifepp::TemMatrix are scalar !
55 */
56
57
58 #ifndef TERM_MATRIX_HPP
59 #define TERM_MATRIX_HPP
60
61 #include "config.h"
62 #include "Term.hpp"
63 #include "TermVector.hpp"
64 #include "SuTermMatrix.hpp"
65 #include "form.h"
66 #include "largeMatrix.h"
67 #include "solvers.h"
68 #include "computation/termUtils.hpp"
69 #include "EigenTerms.hpp"
70
71 namespace xlifepp
72 {
73
74 //@{
75 //! useful aliases for TermMatrix class
76 typedef std::map<uvPair, SuTermMatrix*>::iterator it_mustm;
77 typedef std::map<uvPair, SuTermMatrix*>::const_iterator cit_mustm;
78 //@}
79
80 extern TermMatrix theDefaultTermMatrix;
81 extern const TermVector theDefaultTermVector;
82 extern Preconditioner theDefaultPreconditioner;
83
84 /*!
85 \class TermMatrix
86 end user class handling numerical representation of any
87 bilinear form (BilinearForm). May be multiple unknowns
88 */
89 class TermMatrix : public Term
90 {
91 protected:
92 BilinearForm bilinForm_; //!< bilinear form
93 std::map<uvPair, SuTermMatrix*> suTerms_; //!< list of single unknowns term matrices
94 MatrixEntry* entries_p; //!< pointer to values of TermMatrix (see explanations)
95
96 //to deal with constraints
97 SetOfConstraints* constraints_u_p; //!< constraints data related to conditions to apply to unknowns (col)
98 SetOfConstraints* constraints_v_p; //!< constraints data related to conditions to apply to test functions (row), in general equal to constraints_u_p
99 MatrixEntry* rhs_matrix_p; //!< pointer to a correction matrix when constraints applied
100
101 // for scalar representation
102 MatrixEntry* scalar_entries_p; //!< values of TermMatrix in scalar representation, same as entries_p if scalar unknowns
103 std::vector<DofComponent> cdofs_c; //!< col scalar numbering when scalar_entries_p is allocated
104 std::vector<DofComponent> cdofs_r; //!< row scalar numbering when scalar_entries_p is allocated
105
106 public:
107 TermMatrix(const string_t & na="?"); //!< default constructor
108
109 //old constructors kept for ascending compatibility
110 TermMatrix(const BilinearForm &, const string_t& na, bool noass); //!< constructor with no essential boundary conditions
111 TermMatrix(const BilinearForm &, const EssentialConditions&,
112 const string_t& na, bool noass); //!< constructor with blf and essential conditions on u
113 TermMatrix(const BilinearForm &, const EssentialConditions&, const EssentialConditions&,
114 const string_t& na, bool noass ); //!< constructor with blf and essential conditions on u and v (if different)
115 TermMatrix(const BilinearForm &, const SetOfConstraints&,
116 const string_t& na, bool noass); //!< constructor with blf and explicit constraints data on u
117 TermMatrix(const BilinearForm &, const SetOfConstraints&, const SetOfConstraints&,
118 const string_t& na, bool noass); //!< constructor with blf and explicit constraints data on u and v (if different)
119
120 //constructors with options and no essential boundary conditions
121 TermMatrix(const BilinearForm&, const string_t& na="");
122 TermMatrix(const BilinearForm&, TermOption, const string_t& na="");
123 TermMatrix(const BilinearForm&, TermOption, TermOption, const string_t& na="");
124 TermMatrix(const BilinearForm&, TermOption, TermOption, TermOption, const string_t& na="");
125
126 //constructor with options and one essential boundary condition (same on u and v)
127 TermMatrix(const BilinearForm&, const EssentialConditions&, const string_t& na="");
128 TermMatrix(const BilinearForm&, const EssentialConditions&, TermOption, const string_t& na="");
129 TermMatrix(const BilinearForm&, const EssentialConditions&, TermOption, TermOption, const string_t& na="");
130 TermMatrix(const BilinearForm&, const EssentialConditions&, TermOption, TermOption, TermOption, const string_t& na="");
131
132 //constructor with one essential boundary condition (same on u and v) and explicit reduction method
133 TermMatrix(const BilinearForm&, const EssentialConditions&, const ReductionMethod&, const string_t& na="");
134 TermMatrix(const BilinearForm&, const EssentialConditions&, const ReductionMethod&, TermOption, const string_t& na="");
135 TermMatrix(const BilinearForm&, const EssentialConditions&, const ReductionMethod&, TermOption, TermOption, const string_t& na="");
136 TermMatrix(const BilinearForm&, const EssentialConditions&, const ReductionMethod&, TermOption, TermOption, TermOption, const string_t& na="");
137
138 //constructor with options and two essential boundary conditions (resp. on u and v)
139 TermMatrix(const BilinearForm&, const EssentialConditions&, const EssentialConditions&, const string_t& na="");
140 TermMatrix(const BilinearForm&, const EssentialConditions&, const EssentialConditions&, TermOption, const string_t& na="");
141 TermMatrix(const BilinearForm&, const EssentialConditions&, const EssentialConditions&, TermOption, TermOption, const string_t& na="");
142 TermMatrix(const BilinearForm&, const EssentialConditions&, const EssentialConditions&, TermOption, TermOption, TermOption, const string_t& na="");
143
144 //constructor with options and two essential boundary conditions (resp. on u and v) and explicit reduction method
145 TermMatrix(const BilinearForm&, const EssentialConditions&, const EssentialConditions&, const ReductionMethod&, const string_t& na="");
146 TermMatrix(const BilinearForm&, const EssentialConditions&, const EssentialConditions&, const ReductionMethod&, TermOption, const string_t& na="");
147 TermMatrix(const BilinearForm&, const EssentialConditions&, const EssentialConditions&, const ReductionMethod&, TermOption, TermOption, const string_t& na="");
148 TermMatrix(const BilinearForm&, const EssentialConditions&, const EssentialConditions&, const ReductionMethod&, TermOption, TermOption, TermOption, const string_t& na="");
149
150 void build(const BilinearForm&, const EssentialConditions*, const EssentialConditions*, const ReductionMethod&,
151 const std::vector<TermOption>&, const string_t&); //!< effective constructor from bilinear form, essential conditions and options
152
153 //special constructors
154 TermMatrix(const TermMatrix &, const string_t& na=""); //!< copy constructor (full copy)
155 explicit TermMatrix(const TermVector &, const string_t& na=""); //!< constructor of a diagonal TermMatrix from a TermVector
156 TermMatrix(const TermMatrix&, SpecialMatrix, const string_t& na=""); //!< construct special matrix from a TermMatrix (e.g _IdMatrix)
157 TermMatrix(const SuTermMatrix &, const string_t& na=""); //!< constructor of a TermMatrix from a SuTermMatrix (full copy)
158 TermMatrix(SuTermMatrix *, const string_t& na=""); //!< constructor of a TermMatrix from a SuTermMatrix * (pointer copy)
159 template<typename T>
160 TermMatrix(const Unknown&, const GeomDomain&, const Vector<T> &, const string_t& na=""); //!< special constructor of a diagonal TermMatrix from a vector
161 TermMatrix(const Unknown&, const GeomDomain&, const Unknown&, const GeomDomain&,
162 const OperatorOnKernel&, const string_t& na=""); //!< constructor of matrix opker(xi,yj) (Lagrange interpolation)
163 TermMatrix(const LcTerm<TermMatrix>&,const string_t& na=""); //!<constructor from a linear combination of TermMatrix's
164
165 //destructor, assign and construction tools
166 virtual ~TermMatrix(); //!< destructor
167 TermMatrix& operator=(const LcTerm<TermMatrix>&); //!< assign operator from a linear combination of TermMatrix's
168 TermMatrix& operator=(const TermMatrix&); //!< assign operator (full copy)
169 void initFromBlf(const BilinearForm &, const string_t& na ="", bool noass = false); //!< initializer with blf
170 void copy(const TermMatrix&); //!< copy a TermMatrix in current TermMatrix
171 void clear(); //!< deallocate memory used by a TermMatrix
172 void insert(const SuTermMatrix &); //!< insert a SuTermMatrix with full copy
173 void insert(SuTermMatrix *); //!< insert a SuTermMatrix with pointer copy
174 void initTermVector(TermVector&, ValueType, bool col=true) const; //!< initialize a TermVector consistent with matrix rows or columns
175 template<typename T>
176 void initTermVectorT(TermVector&,const T&, bool col=true) const; //!< initialize a TermVector from value T, consistent with matrix rows or columns
177
178 //accessors and properties
name() const179 string_t name() const { return Term::name(); }
180 void name(const string_t& nam); //!< update the name of a TermMatrix and its SuTermMatrixs
181 void setStorage(StorageType, AccessType); //!< specify the storage of SuTermMatrix's
storageType()182 StorageType storageType() //! returns storage type (cs, skyline, dense)
183 {return computingInfo_.storageType;}
storageAccess()184 AccessType storageAccess() //! returns storage access (row, col, sym, dual)
185 {return computingInfo_.storageAccess;}
isSingleUnknown() const186 bool isSingleUnknown() const //! true if a single unknown matrix
187 {return suTerms_.size()==1;}
188 std::set<const Unknown*> rowUnknowns() const; //!< return set of row Unknowns (say v unknowns)
189 std::set<const Unknown*> colUnknowns() const; //!< return set of col Unknowns (say u unknowns)
190 std::set<const Space*> unknownSpaces() const; //!< list of involved unknown spaces
191 ValueType valueType() const; //!< return value type (_real,_complex)
192 bool isScalar() const; //!< true if all SuTermMatrix's are scalar
hasConstraints() const193 bool hasConstraints() const //! true if one constraints pointer is not null
194 {return (constraints_u_p!=0 || constraints_v_p!=0);}
195 number_t numberOfRows() const; //!< number of rows (counted in scalar dof)
196 number_t numberOfCols() const; //!< number of columns (counted in scalar dof)
197 number_t nnz() const; //!< number of non zero coefficients (counted in scalar)
198 void markAsComputed(bool = true); //!< TermMatrix and SutermMatrix's computed state = true or false
199
200 TermMatrix operator()(const Unknown&,const Unknown&) const; //!< access to single unknowns block matrix as TermMatrix
201 SuTermMatrix& subMatrix(const Unknown*,const Unknown*); //!< access to single unknowns term matrix
202 const SuTermMatrix& subMatrix(const Unknown*,const Unknown*) const; //!< access to single unknowns term matrix (const)
203 SuTermMatrix* subMatrix_p(const Unknown*,const Unknown*); //!< access to single unknowns term matrix pointer
204 const SuTermMatrix* subMatrix_p(const Unknown*,const Unknown*) const; //!< access to single unknowns term matrix pointer (const)
205 TermVector& getRowCol(number_t, AccessType, TermVector&) const; //!< return the i-th row or col as a TermVector (i=1,...)
206 TermVector row(number_t) const; //!< return the i-th row as a TermVector (i=1,...)
207 TermVector column(number_t) const; //!< return the j-th col as a TermVector (j=1,...)
208 template <typename T>
209 LargeMatrix<T>& getLargeMatrix(StrucType =_undefStrucType) const; //!< get internal matrix (restrictive), advanced usage
210 template <typename T>
211 HMatrix<T,FeDof>& getHMatrix(StrucType =_undefStrucType) const; //!< get internal HMatrix if exists, advanced usage
212
begin() const213 cit_mustm begin() const //! iterator to the first element of the SuTermMatrix map (const)
214 {return suTerms_.begin();}
begin()215 it_mustm begin() //! iterator to the first element of the SuTermMatrix map
216 {return suTerms_.begin();}
end() const217 cit_mustm end() const //! iterator to the last element of the SuTermMatrix map (const)
218 {return suTerms_.end();}
end()219 it_mustm end() //! iterator to the last element of the SuTermMatrix map
220 {return suTerms_.end();}
firstSut() const221 SuTermMatrix* firstSut() const //! first SuTermMatrix as pointer, 0 if void TermMatrix (const)
222 {
223 if(suTerms_.size()>0) return suTerms_.begin()->second;
224 else return 0;
225 }
firstSut()226 SuTermMatrix* firstSut() //! first SuTermMatrix as pointer, 0 if void TermMatrix
227 {
228 if(suTerms_.size()>0) return suTerms_.begin()->second;
229 else return 0;
230 }
suTerms()231 std::map<uvPair, SuTermMatrix*>& suTerms() //! access to suterms
232 {return suTerms_;}
nbTerms() const233 number_t nbTerms() const //! number of suTerms
234 {return suTerms_.size();}
entries()235 MatrixEntry*& entries() //! access to entries pointer (r/w)
236 {return entries_p;}
entries() const237 const MatrixEntry* entries() const //! access to entries pointer (r)
238 {return entries_p;}
scalar_entries()239 MatrixEntry*& scalar_entries() //! access to scalar entries pointer (r/w)
240 {return scalar_entries_p;}
scalar_entries() const241 const MatrixEntry* scalar_entries() const //! access to scalar entries pointer (r)
242 {return scalar_entries_p;}
243 MatrixEntry* actual_entries() const; //!< return the actual pointer to entries (priority to scalar entry)
rhs_matrix()244 MatrixEntry*& rhs_matrix() //! access to rhs matrix pointer (r/w)
245 {return rhs_matrix_p;}
rhs_matrix() const246 const MatrixEntry* rhs_matrix() const //! access to rhs matrix pointer (r)
247 {return rhs_matrix_p;}
cdofsr() const248 const std::vector<DofComponent>& cdofsr() const //! access to cdofs_r vector (r)
249 {return cdofs_r;}
cdofsc() const250 const std::vector<DofComponent>& cdofsc() const //! access to cdofs_c vector (r)
251 {return cdofs_c;}
cdofsr()252 std::vector<DofComponent>& cdofsr() //! access to cdofs_r vector (r/w)
253 {return cdofs_r;}
cdofsc()254 std::vector<DofComponent>& cdofsc() //! access to cdofs_c vector (r/w)
255 {return cdofs_c;}
256
257 SymType symmetry() const; //!< return global symmetry property
258 FactorizationType factorization() const; //!< return factorization type if factorized
259
260 //computation functions and related tools
261 void compute(); //!< compute TermMatrix
262 void compute(const LcTerm<TermMatrix>&,const string_t& na=""); //!< compute TermMatrix from a linear combination of TermMatrix's
263
264 //storage management tools
265 std::pair<StorageType,AccessType> findGlobalStorageType() const; //!< suggest storage type for global representation
266 void toScalar(bool keepEntries=false); //!< create scalar representation of TermMatrix
267 void toGlobal(StorageType, AccessType, SymType=_noSymmetry,
268 bool keepSuTerms=false); //!< create global scalar representation of TermMatrix
269 void toSkyline(); //!< convert to skyline storage (works only for cs storage single unknown matrix)
270 void mergeNumbering(std::map<const Unknown*, std::list<SuTermMatrix* > >&, std::map<SuTermMatrix*, std::vector<number_t > >&,
271 std::vector<DofComponent>&, AccessType); //!< merging tool used by toGlobal()
272 void mergeNumberingFast(std::map<const Unknown*, std::list<SuTermMatrix*> >&, std::map<SuTermMatrix*, std::vector<number_t > >&,
273 std::vector<DofComponent>&, AccessType); //!< merging tool used by toGlobal(), fast version
274
275 void addIndices(std::vector<std::set<number_t> >&, MatrixStorage*,
276 const std::vector<number_t>&, const std::vector<number_t>&); //!< merging tool used by toGlobal()
277 void mergeBlocks(); //!< merge blocks referring to same vector unknowns
278 MatrixEntry* matrixData(); //!< turns to scalar representation and returns pointer to data
279
280 //essential conditions tools
281 void pseudoReduction(); //!< reduction of essential condition using pseudo reduction method
282
283 //scalar type moving
284 TermMatrix& toReal(); //!< go to real representation getting the real part when complex
285 TermMatrix& toComplex(); //!< go to complex representation
286 TermMatrix& toImag(); //!< go to real representation getting the imag part when complex
287 TermMatrix& toConj(); //!< change to conj(U)
288 TermMatrix& roundToZero(real_t aszero=10*theEpsilon); //!< round to zero all coefficients close to 0 (|.| < aszero)
289 friend TermMatrix real(const TermMatrix&);
290 friend TermMatrix imag(const TermMatrix&);
291 friend TermMatrix toComplex(const TermMatrix&);
292 friend TermMatrix roundToZero(const TermMatrix&, real_t aszero);
293
294 real_t norm2() const; //!< Frobenius norm : sqrt(sum|a_ij|^2)
295 real_t norminfty() const; //!< infinite norm : sup|a_ij|
296
297 //in/out utilities
298 void print(std::ostream&) const; //!< print TermMatrix
print(PrintStream & os) const299 void print(PrintStream& os) const {print(os.currentStream());}
300 void printSummary(std::ostream&) const; //!< print summary TermMatrix
printSummary(PrintStream & os) const301 void printSummary(PrintStream& os) const {printSummary(os.currentStream());}
302 void viewStorage(std::ostream&) const; //!< print storage on stream
viewStorage(PrintStream & os) const303 void viewStorage(PrintStream& os) const {viewStorage(os.currentStream());}
304 void saveToFile(const string_t&, StorageType, bool=false) const; //!< save TermMatrix to file (dense or Matlab format)
305
306 //various algebraic operation
307 template<typename T> TermMatrix& operator*=(const T&); //!< TermMatrix *= t
308 template<typename T> TermMatrix& operator/=(const T&); //!< TermMatrix /= t
309 TermMatrix& operator+=(const TermMatrix&); //!< TermMatrix += TermMatrix
310 TermMatrix& operator-=(const TermMatrix&); //!< TermMatrix -= TermMatrix
311 friend TermVector& multMatrixVector(const TermMatrix&,const TermVector&, TermVector&);
312 friend TermVector operator*(const TermMatrix&,const TermVector&);
313 friend TermVector& multVectorMatrix(const TermVector&, const TermMatrix&, TermVector&);
314 friend TermVector& multVectorMatrix(const TermMatrix&, const TermVector&, TermVector&);
315 friend TermVector operator*(const TermVector&, const TermMatrix&);
316 friend TermMatrix operator*(const TermMatrix&, const TermMatrix&);
317
318 //solvers
319 friend TermVector prepareLinearSystem(TermMatrix&, TermVector&, MatrixEntry*&, VectorEntry*&,
320 StorageType, AccessType, bool);
321 friend void updateRhs(TermMatrix&, TermVector&);
322 void ldltSolve(const TermVector&, TermVector&); //!< solve linear system using LDLt factorization
323 void ldlstarSolve(const TermVector&,TermVector&); //!< solve linear system using LDL* factorization
324 void luSolve(const TermVector&,TermVector&); //!< solve linear system using LU factorization
325 void umfpackSolve(const TermVector&,TermVector&); //!< solve linear system using umfpack
326 void sorSolve(const TermVector& V, TermVector& R, const real_t w, SorSolverType sType ); //!< solve linear system using SOR
327 void sorUpperSolve(const TermVector& V, TermVector& R, const real_t w ); //!< solve linear system using SOR for upper triangular matrices
328 void sorDiagonalMatrixVector(const TermVector& V, TermVector& R, const real_t w ); //!< solve linear system using SOR for block diagonal matrices
329 void sorLowerSolve(const TermVector& V, TermVector& R, const real_t w ); //!< solve linear system using SOR for lower triangular matrices
330 void sorDiagonalSolve(const TermVector& V, TermVector& R, const real_t w ); //!< solve linear system using SOR for diagonal matrices
331 };
332
333 TermMatrix real(const TermMatrix&); //!< return real part as a real TermMatrix
334 TermMatrix imag(const TermMatrix&); //!< return imag part as a real TermMatrix
335 TermMatrix toComplex(const TermMatrix&); //!< return the complex representation of the given TermMatrix
336 TermMatrix roundToZero(const TermMatrix&, real_t aszero); //!< return the 0 rounded TermMatrix
337
338 TermVector& multMatrixVector(const TermMatrix&,const TermVector&, TermVector&); //!< product TermMatrix * TermVector
339 TermVector operator*(const TermMatrix&,const TermVector&); //!< product TermMatrix * TermVector
340 TermVector& multVectorMatrix(const TermVector&, const TermMatrix&, TermVector&); //!< product TermVector * TermMatrix
341 TermVector& multVectorMatrix(const TermMatrix&, const TermVector&, TermVector&); //!< product TermVector * TermMatrix
342 TermVector operator*(const TermVector&, const TermMatrix&); //!< product TermVector * TermMatrix
343 TermMatrix operator*(const TermMatrix&, const TermMatrix&); //!< product of TermMatrix
344
345 //! prepare linear system AX=B (internal tool)
346 TermVector prepareLinearSystem(TermMatrix&, TermVector&, MatrixEntry*&, VectorEntry*&, StorageType, AccessType, bool);
347 void updateRhs(TermMatrix&, TermVector&); //!< prepare rhs B of linear system AX=B (internal tool)
348
349 // ------------------------------------------------------------------------------------------------------------------------
350 // External functions declaration
351 // ------------------------------------------------------------------------------------------------------------------------
352 #include "decLinSys.hpp"
353 #include "decEigen.hpp"
354
norm2(const TermMatrix & A)355 inline real_t norm2(const TermMatrix& A)
356 {return A.norm2();}
357
norminfty(const TermMatrix & A)358 inline real_t norminfty(const TermMatrix& A)
359 {return A.norminfty();}
360
361 //in/out utils
saveToFile(const TermMatrix & A,const string_t & fn,StorageType st,bool enc=false)362 inline void saveToFile(const TermMatrix& A, const string_t& fn, StorageType st, bool enc = false) //! save TermMatrix to file (dense or Matlab format)
363 {A.saveToFile(fn,st,enc);}
364
365 //integral representation
366 TermMatrix integralRepresentation(const Unknown&, const GeomDomain&, const LinearForm&, string_t nam="IR"); //!< integral representation
367 TermMatrix integralRepresentation(const GeomDomain&, const LinearForm&, string_t nam="IR"); //!< integral representation, no unknown
368 TermMatrix integralRepresentation(const std::vector<Point>&, const LinearForm&, string_t nam="IR"); //!< integral representation, no unknown, no domain
369
370 //----------------------------------------------------------------------------------------------------------
371 // template functions
372 //----------------------------------------------------------------------------------------------------------
373 //! special constructor of a diagonal TermMatrix from a vector
374 template<typename T>
TermMatrix(const Unknown & u,const GeomDomain & dom,const Vector<T> & v,const string_t & na)375 TermMatrix::TermMatrix(const Unknown& u, const GeomDomain& dom, const Vector<T> & v, const string_t& na)
376 {
377 trace_p->push("TermMatrix::TermMatrix(Vector<T>)");
378 entries_p = 0;
379 scalar_entries_p=0;
380 rhs_matrix_p=0;
381 constraints_u_p = 0;
382 constraints_v_p = 0;
383 termType_ = _termMatrix;
384 computingInfo_.noAssembly = false;
385 name_ = na;
386 TermVector tv(u,dom,v);
387 for(it_mustv it = tv.begin(); it != tv.end(); it++)
388 suTerms_[uvPair(&u,&u)]=new SuTermMatrix(*it->second, na+"_"+u.name());
389 computingInfo_.isComputed = true;
390 trace_p->pop();
391 }
392
393 /*
394 multiple algebraic operations on TermMatrix
395 produce a LcTerm object (linear combination of terms)
396 the computation is done by constructor from LcTerm or assign operation of a LcTerm
397 */
398 const TermMatrix& operator+(const TermMatrix&); //!< unary operator+
399 LcTerm<TermMatrix> operator-(const TermMatrix&); //!< unary operator- (returns a LcTerm)
400 LcTerm<TermMatrix> operator+(const TermMatrix&, const TermMatrix&); //!< addition of TermMatrix
401 LcTerm<TermMatrix> operator-(const TermMatrix&, const TermMatrix&); //!< substraction of TermMatrix of a TermMatrix
402 LcTerm<TermMatrix> operator+(const LcTerm<TermMatrix>&, const TermMatrix&); //!< addition of a TermMatrix to a LcTerm
403 LcTerm<TermMatrix> operator-(const LcTerm<TermMatrix>&, const TermMatrix&); //!< substraction of a TermMatrix of a LcTerm
404 LcTerm<TermMatrix> operator+(const TermMatrix&, const LcTerm<TermMatrix>&); //!< addition of a TermMatrix to a LcTerm
405 LcTerm<TermMatrix> operator-(const TermMatrix&, const LcTerm<TermMatrix>&); //!< substraction of a LcTerm from a TermMatrix
406 TermVector operator*(const LcTerm<TermMatrix>&, const TermVector&); //!< LcTerm * TermVector
407
408 //! product of TermMatrix by a real or a complex (template T)
409 template<typename T>
operator *(const TermMatrix & tv,const T & t)410 LcTerm<TermMatrix> operator*(const TermMatrix& tv, const T& t)
411 {
412 return LcTerm<TermMatrix>(&tv,t);
413 }
414 //! product of TermMatrix by a real or a complex (template T)
415 template<typename T>
operator *(const T & t,const TermMatrix & tv)416 LcTerm<TermMatrix> operator*(const T& t, const TermMatrix& tv)
417 {
418 return LcTerm<TermMatrix>(&tv,t);
419 }
420 //! division of TermMatrix by a real or a complex (template T)
421 template<typename T>
operator /(const TermMatrix & tv,const T & t)422 LcTerm<TermMatrix> operator/(const TermMatrix& tv, const T& t)
423 {
424 return LcTerm<TermMatrix>(&tv,1./t);
425 }
426
427 //! product assignment by a real or a complex (template T)
428 template<typename T>
operator *=(const T & t)429 TermMatrix& TermMatrix::operator*=(const T& t)
430 {
431 for(it_mustm it = suTerms_.begin(); it != suTerms_.end(); it++) *it->second *= t;
432 bilinForm_*=t;
433 return *this;
434 }
435
436 //! division assignment by a real or a complex (template T)
437 template<typename T>
operator /=(const T & t)438 TermMatrix& TermMatrix::operator/=(const T& t)
439 {
440 if(t == 0) error("divBy0");
441 for(it_mustm it = suTerms_.begin(); it != suTerms_.end(); it++) *it->second /= t;
442 bilinForm_/=t;
443 return *this;
444 }
445
446 //! product of TermMatrix by a real or a complex (template T)
447 template<typename T>
operator *(const TermMatrix & tM,const T & t)448 TermMatrix operator*(const TermMatrix & tM, const T& t)
449 {
450 TermMatrix R(tM);
451 return R*=t;
452 }
453
454 //! initialize a TermVector consistent with matrix column unknown (col=true), with matrix row unknown (col=false)
455 template<typename T>
initTermVectorT(TermVector & tv,const T & v,bool col) const456 void TermMatrix::initTermVectorT(TermVector& tv, const T& v, bool col) const
457 {
458 string_t na = tv.name(); //preserving only name
459 tv.clear();
460 tv.name() = na;
461 //make list of col unknowns
462 std::map<const Unknown*, SuTermMatrix*> terms;
463 for(cit_mustm itA = begin(); itA != end(); itA++)
464 {
465 const Unknown* u;
466 if(col) u = itA->first.first;
467 else u = itA->first.second;
468 SuTermMatrix* Auv = itA->second;
469 if(terms.find(u) == terms.end()) terms[u] = Auv;
470 }
471 std::map<const Unknown*, SuTermMatrix*>::iterator it;
472 for(it = terms.begin(); it != terms.end(); it++)
473 {
474 const Unknown* u = it->first;
475 SuTermMatrix* sut = it->second;
476 Space* sp = sut->space_up();
477 if(!col) sp = sut->space_vp();
478 SuTermVector* sutX = new SuTermVector("", u, sp, v);
479 tv.insert(u, sutX);
480 }
481 }
482
483 /*! get internal matrix (restrictive), i.e LargeMatrix object if available
484 available only for single unknowns TermMatrix and multiple unknowns TermMatrix having a global representation
485 if scalar and vector representations both exist, the scalar representation is prior except if StrucType argument is specified
486 user has to specify the LargeMatrix type :
487 TermMatrix M(intg(omega,u*v));
488 LargeMatrix& Alm=M.getMatrix<Real>();
489 It returns a reference to the LargeMatrix
490 Note : designed for advanced usage
491 */
492 template <typename T>
getLargeMatrix(StrucType str) const493 LargeMatrix<T>& TermMatrix::getLargeMatrix(StrucType str) const
494 {
495 number_t nbsut=suTerms_.size();
496 if (nbsut==1) return suTerms_.begin()->second->getLargeMatrix<T>(str);
497 if (scalar_entries_p!=0 && str!=_matrix) return scalar_entries_p->getLargeMatrix<T>();
498 if (entries_p!=0) return entries_p->getLargeMatrix<T>();
499 where("TermMatrix::getLargeMatrix(StrucType)");
500 error("term_no_entries");
501 return *(new LargeMatrix<T>()); //fake return
502 }
503
504 /*! get internal H-matrix (restrictive), i.e HMatrix object if available
505 available only for single unknown TermMatrix
506 if scalar and vector representations both exist, the scalar representation is prior except if StrucType argument is specified
507 It returns a reference to a HMatrix
508 Note : designed for advanced usage
509 */
510 template <typename T>
getHMatrix(StrucType str) const511 HMatrix<T,FeDof>& TermMatrix::getHMatrix(StrucType str) const
512 {
513 number_t nbsut=suTerms_.size();
514 if (nbsut==1) return suTerms_.begin()->second->getHMatrix<T>(str);
515 where("TermMatrix::getHMatrix(StrucType)");
516 error("term_no_entries");
517 return *(new HMatrix<T,FeDof>()); //fake return
518 }
519
520 /* ================================================================================================
521 computation algorithm for bilinear form as linear form
522 (placed here because it uses TermMatrix's machinery)
523 ================================================================================================ */
524 /*!
525 Computation of a linear form defined from a bilinear form and a TermVector to apply to
526 l(v) = a(U,v) or l(u)=a(u,V) where U,V are some term vectors
527 \param lf : a pair of BasicLinearForm and coefficient
528 \param val : vector of values to fill in
529 \param vt : type of value type
530
531 example : intg_gamma grad(U)|grad(v) dy U a termvector
532 int_sigma intg_gamma (U*G)*v U a termvector
533
534 The algorithm is the most lazy one :
535 - go to the TermMatrix machinery to compute matrix A from a(u,v)
536 - do the product A*U or V*A
537 - add the product result to vector val
538 */
539
540 template<typename T, typename K>
computeBilAsLin(const std::pair<BasicLinearForm *,complex_t> & lf,Vector<T> & val,K & vt)541 void computeBilAsLin(const std::pair<BasicLinearForm*, complex_t>& lf, Vector<T>& val, K& vt)
542 {
543 if(lf.first->type()!=_bilinearAsLinear)
544 {
545 where("SuTermVector::computeBilAsLin(...)");
546 error("lform_not_handled", words("form type", lf.first->type()));
547 }
548 BilinearFormAsLinearForm* blfaslf = lf.first->asBilinearAsLinearForm();
549 K coef = complexToT<K>(lf.second);
550 //create TermMatrix
551 SuBilinearForm sublf(std::vector<blfPair>(1, blfPair(blfaslf->blf(), lf.second)));
552 BilinearForm blf(sublf);
553 TermMatrix S(blf);
554 //do product
555 TermVector R;
556 if(blfaslf->asUnknown_) R=S* *blfaslf->termVector();
557 else R=*blfaslf->termVector()*S;
558 //add R to vector val
559 Vector<T> r;
560 R.asVector(r);
561 typename Vector<T>::iterator itv=val.begin(), itr=r.begin();
562 for(;itr!=r.end(); ++itr, ++itv) *itv+=*itr;
563 }
564
565 } //end of namespace xlifepp
566
567 #endif /* TERM_MATRIX_HPP */
568