G+Smo  25.01.0
Geometry + Simulation Modules
 
Loading...
Searching...
No Matches
gsSparseSystem.h
Go to the documentation of this file.
1
15#pragma once
16
17#include <gsCore/gsStdVectorRef.h>
18
19namespace gismo
20{
21
28template<typename T, bool symm>
29class gsSparseSystem
30{
31protected:
32 typedef gsStdVectorRef<gsDofMapper> DofMapperRef;
33
34 typedef std::vector<gsDofMapper> DofMappers;
35
36 typedef typename gsSparseMatrix<T>::BlockView matBlockView;
37
38 typedef typename gsMatrix<T>::BlockView rhsBlockView;
39
40protected:
41
43 gsSparseMatrix<T> m_matrix;
44
46 gsMatrix<T> m_rhs;
47
48
49 // -- Structure
50
56 DofMappers m_mappers;
57
60 gsVector<index_t> m_row;
61
64 gsVector<index_t> m_col;
65
68 gsVector<index_t> m_rstr;
69
72 gsVector<index_t> m_cstr;
73
77 gsVector<index_t> m_cvar;
78
79 gsVector<index_t> m_dims;
80
81public:
82
83 gsSparseSystem()
84 { }
85
91 gsSparseSystem(gsDofMapper & mapper)
92 : m_mappers(1),
93 m_row (1),
94 m_col (1),
95 m_rstr (1),
96 m_cstr (1),
97 m_cvar (1),
98 m_dims (1)
99 {
100 m_row [0] = m_col [0] =
101 m_rstr[0] = m_cstr[0] =
102 m_cvar[0] = 0;
103
104 m_dims[0] = 1;
105
106 m_mappers.front().swap(mapper);
107
108 m_matrix.resize( m_mappers.front().freeSize() ,
109 m_mappers.front().freeSize() );
110 }
111
112
139 gsSparseSystem(DofMappers & mappers,
140 const gsVector<index_t> & dims)
141 : m_row(dims.sum()),
142 m_col(dims.sum()),
143 m_rstr(dims.sum()),
144 m_cstr(dims.sum()),
145 m_dims(dims.cast<index_t>())
146 {
147 const index_t d = dims.size();
148 const index_t s = dims.sum();
149 const index_t ms = mappers.size();
150
151 GISMO_ASSERT(ms==s ||ms==2*s, "Connot deduce block structure");
152
153 m_mappers.swap(mappers);
154
155 //Calculate the map for blocks to mappers
156 index_t k=0;
157 for(index_t i=0;i<d;++i)
158 for(index_t j=0; j<dims[i];++j)
159 {
160 m_row[k]=i;
161 ++k;
162 }
163 //At first glance, row blocks and column blocks have the same mapper
164 m_col = m_row;
165
166 if (ms == 1 )
167 {
168 m_row.setZero();
169 m_col.setZero();
170 m_cvar.setZero(1);
171 }
172 else if ( ms == 2*s )
173 m_col.array() += s; //mappers s+1... 2s are then the column mappers
174
175 //assumes that the mappers are ordered as the bases in gsAssembler and are starting from m_bases[0]!
176 m_cvar = m_row.cast<index_t>();
177
178 //Fill the strides
179 m_rstr[0] = m_cstr[0] = 0;
180 for (index_t r = 1; r < d; ++r) // for all row-blocks
181 m_rstr[r] = m_rstr[r-1] + m_mappers[m_row[r-1]].freeSize();
182 for (index_t c = 1; c < d; ++c) // for all col-blocks
183 m_cstr[c] = m_cstr[c-1] + m_mappers[m_col[c-1]].freeSize();
184
185 m_matrix.resize( m_rstr.at(d-1) + m_mappers[m_row[d-1]].freeSize() ,
186 m_cstr.at(d-1) + m_mappers[m_col[d-1]].freeSize() );
187
188 }
189
205 gsSparseSystem(DofMappers & mappers,
206 const index_t rows,
207 const index_t cols)
208 : m_row (gsVector<index_t>::LinSpaced(rows,0,rows-1)),
209 m_col (gsVector<index_t>::LinSpaced(cols,0,cols-1)),
210 m_rstr(rows),
211 m_cstr(cols),
212 m_dims(cols)
213 {
214 GISMO_ASSERT( rows > 0 && cols > 0, "Block dimensions must be positive");
215
216 m_mappers.swap(mappers);
217
218 if ( static_cast<index_t>(m_mappers.size()) == rows + cols )
219 {
220 m_col.array() += rows;
221 m_cvar = m_row.cast<index_t>();//assume 1-1 bases
222 }
223 else if ( static_cast<index_t>(m_mappers.size()) == rows )
224 {
225 GISMO_ENSURE( rows == cols, "Dof Mapper vector does not match block dimensions");
226 m_cvar = m_row.cast<index_t>();//assume 1-1 bases
227 }
228 else if ( m_mappers.size() == 1 )
229 {
230 m_row.setZero();
231 m_col.setZero();
232 m_cvar.setZero(1);
233 }
234 else
235 {
236 GISMO_ERROR("Cannot deduce block structure.");
237 }
238
239 m_dims.setOnes();
240
241 m_rstr[0] = m_cstr[0] = 0;
242 for (index_t r = 1; r < rows; ++r) // for all row-blocks
243 m_rstr[r] = m_rstr[r-1] + m_mappers[m_row[r-1]].freeSize();
244 for (index_t c = 1; c < cols; ++c) // for all col-blocks
245 m_cstr[c] = m_cstr[c-1] + m_mappers[m_col[c-1]].freeSize();
246
247 m_matrix.resize( m_rstr.at(rows-1) + m_mappers[m_row[rows-1]].freeSize() ,
248 m_cstr.at(cols-1) + m_mappers[m_col[cols-1]].freeSize() );
249 }
250
275 gsSparseSystem(DofMappers & mappers,
276 const gsVector<index_t> & rowInd,
277 const gsVector<index_t> & colInd,
278 const gsVector<index_t> & colvar)
279 : m_row (rowInd),
280 m_col (colInd),
281 m_rstr((index_t)rowInd.size()),
282 m_cstr((index_t)colInd.size()),
283 m_dims(colInd.size())
284 // ,m_cvar(colvar) //<< Bug
285 {
286 m_dims.setOnes();
287 m_cvar = colvar;
288 const index_t rows = m_row.size();
289 const index_t cols = m_col.size();
290 GISMO_ASSERT( rows > 0 && cols > 0,
291 "Block dimensions must be positive");
292
293 m_mappers.swap(mappers);
294
295 m_rstr[0] = m_cstr[0] = 0;
296 for (index_t r = 1; r < rows; ++r) // for all row-blocks
297 m_rstr[r] = m_rstr[r-1] + m_mappers[m_row[r-1]].freeSize();
298 for (index_t c = 1; c < cols; ++c) // for all col-blocks
299 m_cstr[c] = m_cstr[c-1] + m_mappers[m_col[c-1]].freeSize();
300
301 m_matrix.resize( m_rstr.at(rows-1) + m_mappers[m_row[rows-1]].freeSize() ,
302 m_cstr.at(cols-1) + m_mappers[m_col[cols-1]].freeSize() );
303 }
304
309 void swap(gsSparseSystem & other)
310 {
311 m_matrix .swap(other.m_matrix );
312 m_rhs .swap(other.m_rhs );
313 m_mappers.swap(other.m_mappers);
314 m_row .swap(other.m_row );
315 m_col .swap(other.m_col );
316 m_rstr .swap(other.m_rstr );
317 m_cstr .swap(other.m_cstr );
318 m_cvar .swap(other.m_cvar );
319 m_dims .swap(other.m_dims );
320 }
321
327 void reserve(const index_t nz, const index_t numRhs)
328 {
329 GISMO_ASSERT( 0 != m_mappers.size(), "Sparse system was not initialized");
330 if ( 0 != m_matrix.cols() )
331 {
332 m_matrix.reservePerColumn(nz);
333 if ( 0 != numRhs )
334 m_rhs.setZero(m_matrix.cols(), numRhs);
335 }
336 }
337
350 void reserve(const gsMultiBasis<T> & mb, const gsOptionList & opt,
351 const index_t numRhs)
352 {
353 reserve(numColNz(mb,opt), numRhs);
354 }
355
359 index_t numColNz(const gsMultiBasis<T> & mb, const gsOptionList & opt) const
360 {
361 // Pick up values from options
362 const T bdA = opt.getReal("bdA");
363 const index_t bdB = opt.getInt("bdB");
364 const T bdO = opt.getReal("bdO");
365 const gsBasis<T> & b = mb[0];
366 T nz = 1;
367 for (short_t i = 0; i != b.dim(); ++i)
368 nz *= bdA * cast<short_t,T>(b.degree(i)) + cast<index_t,T>(bdB);
369 return cast<T,short_t>(nz*(1.0+bdO));
370 }
371
373 void setZero()
374 {
375 m_matrix.setZero();
376 m_rhs .setZero();
377 }
378
380 index_t cols() const { return m_matrix.cols(); }
381
383 index_t rows() const { return m_matrix.rows(); }
384
386 index_t rhsRows() const { return m_rhs.rows(); }
387
389 index_t rhsCols() const { return m_rhs.cols(); }
390
391public: /* Accessors */
392
394 const gsSparseMatrix<T> & matrix() const
395 { return m_matrix; }
396
398 gsSparseMatrix<T> & matrix()
399 { return m_matrix; }
400
402 const gsMatrix<T> & rhs() const
403 { return m_rhs; }
404
406 gsMatrix<T> & rhs()
407 { return m_rhs; }
408
410 matBlockView blockView()
411 {
412 gsVector<index_t> rowSizes(m_row.size()), colSizes(m_col.size());
413
414 for (index_t r = 0; r != rowSizes.size(); ++r) // for all row-blocks
415 rowSizes[r] = m_mappers[m_row[r]].freeSize();
416
417 for (index_t c = 0; c != colSizes.size(); ++c) // for all col-blocks
418 colSizes[c] = m_mappers[m_col[c]].freeSize();
419
420 return m_matrix.blockView(rowSizes,colSizes);
421 }
422
438 matBlockView blockView(size_t numRowBlocksNew, size_t numColBlocksNew, const gsVector<index_t>& rowBlocksNew, const gsVector<index_t>& colBlocksNew)
439 {
440 gsVector<index_t> rowSizes(numRowBlocksNew), colSizes(numColBlocksNew);
441 rowSizes.setZero();
442 colSizes.setZero();
443
444 for (index_t r = 0; r != m_row.size(); ++r) // for all row-blocks
445 rowSizes[rowBlocksNew[r]] += m_mappers[m_row[r]].freeSize();
446
447 for (index_t c = 0; c != m_col.size(); ++c) // for all col-blocks
448 colSizes[colBlocksNew[c]] += m_mappers[m_col[c]].freeSize();
449
450 return m_matrix.blockView(rowSizes, colSizes);
451 }
452
459 rhsBlockView blockViewRhs(size_t numRowBlocksNew, const gsVector<index_t>& rowBlocksNew)
460 {
461 gsVector<index_t> rowSizes(numRowBlocksNew), colSizes(1);
462 rowSizes.setZero();
463 colSizes[0]=1;
464
465 for (index_t r = 0; r != m_row.size(); ++r) // for all row-blocks
466 rowSizes[rowBlocksNew[r]] += m_mappers[m_row[r]].freeSize();
467
468 return m_rhs.blockView(rowSizes, colSizes);
469 }
470
472 index_t numColBlocks() const {return m_col.size();}
473
475 index_t numRowBlocks() const {return m_row.size();}
476
482 const gsDofMapper & rowMapper(const index_t r) const
483 { return m_mappers[m_row[r]]; }
484
490 gsDofMapper & rowMapper(const index_t r)
491 { return m_mappers[m_row[r]]; }
492
498 const gsDofMapper & colMapper(const index_t c) const
499 {
500 return m_mappers[m_col[c]];
501 }
502
508 gsDofMapper & colMapper(const index_t c)
509 { return m_mappers[m_col[c]]; }
510
516 index_t colBasis(const index_t c) const // better name ?
517 { return m_cvar[c]; }
518
520 index_t unkSize(const index_t unk) const
521 {return m_dims[unk];}
522
524 index_t numUnknowns() const {return m_dims.size(); }
525
528 const DofMappers & dofMappers() const
529 { return m_mappers; }
530
532 bool initialized() const
533 {
534 return 0 != m_col.size();
535 }
536
538 bool symmetry() const
539 {
540 return symm;
541 }
542
543 /*
544 gsRefVector<gsDofMapper> colMappers()
545 {
546 return gsRefVector<gsDofMapper>(m_mappers, m_col);
547 }
548
549 gsRefVector<gsDofMapper> rowMappers()
550 {
551 return gsRefVector<gsDofMapper>(m_mappers, m_row);
552 }
553 */
554
555
556public: /* mapping patch-local to global indices */
557
567 void mapRowIndices(const gsMatrix<index_t> & actives,
568 const index_t patchIndex,
569 gsMatrix<index_t> & result,
570 const index_t r = 0) const
571 {
572 m_mappers[m_row.at(r)].localToGlobal(actives, patchIndex, result);
573 }
574
584 void mapColIndices(const gsMatrix<index_t> & actives,
585 const index_t patchIndex,
586 gsMatrix<index_t> & result,
587 const index_t c = 0) const
588 {
589 m_mappers[m_col.at(c)].localToGlobal(actives, patchIndex, result);
590 }
591
600 void mapToGlobalRowIndex(const index_t active,
601 const index_t patchIndex,
602 index_t & result,
603 const index_t r = 0) const
604 {
605 result = m_mappers[m_row.at(r)].index(active, patchIndex)+m_rstr[r];
606 }
607
616 void mapToGlobalColIndex(const index_t active,
617 const index_t patchIndex,
618 index_t & result,
619 const index_t c = 0) const
620 {
621 result = m_mappers[m_col.at(c)].index(active, patchIndex)+m_cstr[c];
622 }
623
624public: /* Add local contributions to system matrix */
625
638 void pushToMatrix(const gsMatrix<T> & localMat,
639 const gsMatrix<index_t> & actives,
640 const gsMatrix<T> & eliminatedDofs,
641 const size_t r = 0, const size_t c = 0)
642 {
643 const index_t numActive = actives.rows();
644 const gsDofMapper & rowMap = m_mappers[m_row.at(r)];
645 GISMO_ASSERT( &rowMap == &m_mappers[m_col.at(c)], "Error");
646
647 for (index_t i = 0; i != numActive; ++i)
648 {
649 const index_t ii = m_rstr.at(r) + actives.at(i); // N_i
650
651 if ( rowMap.is_free_index(actives.at(i)) )
652 {
653 for (index_t j = 0; j != numActive; ++j)
654 {
655 const index_t jj = m_cstr.at(c) + actives.at(j); // N_j
656
657 if ( rowMap.is_free_index( actives.at(j)) )
658 {
659 // If matrix is symmetric, we store only lower
660 // triangular part
661 if ( (!symm) || jj <= ii )
662 m_matrix.coeffRef(ii, jj) += localMat(i, j);
663 }
664 else if(0!=eliminatedDofs.size())
665 {
666 m_rhs.row(ii).noalias() -= localMat(i, j) *
667 eliminatedDofs.row( rowMap.global_to_bindex(actives.at(j)) );
668 }
669 }
670 }
671 }
672 }
673
689 void pushToMatrix(const gsMatrix<T> & localMat,
690 const gsMatrix<index_t> & actives_i,
691 const gsMatrix<index_t> & actives_j,
692 const gsMatrix<T> & eliminatedDofs_j,
693 const size_t r = 0, const size_t c = 0)
694 {
695 const index_t numActive_i = actives_i.rows();
696 const index_t numActive_j = actives_j.rows();
697
698 const gsDofMapper & rowMap = m_mappers[m_row.at(r)];
699 const gsDofMapper & colMap = m_mappers[m_col.at(c)];
700
701 for (index_t i = 0; i != numActive_i; ++i)
702 {
703 const index_t ii = m_rstr.at(r) + actives_i.at(i); // N_i
704
705 if ( rowMap.is_free_index(actives_i.at(i)) )
706 {
707 for (index_t j = 0; j != numActive_j; ++j)
708 {
709 const index_t jj = m_cstr.at(c) + actives_j.at(j); // N_j
710
711 if ( colMap.is_free_index(actives_j.at(j)) )
712 {
713 // If matrix is symmetric, we store only lower
714 // triangular part
715 if ( (!symm) || jj <= ii )
716 m_matrix.coeffRef(ii, jj) += localMat(i, j);
717 }
718 else
719 {
720 m_rhs.row(ii).noalias() -= localMat(i, j) *
721 eliminatedDofs_j.row( colMap.global_to_bindex(actives_j.at(j)) );
722 }
723 }
724 }
725 }
726 }
727
728
740 void pushToMatrixAllFree(const gsMatrix<T> & localMat,
741 const gsMatrix<index_t> & actives,
742 const size_t r = 0, const size_t c = 0)
743 {
744 const index_t numActive = actives.rows();
745 GISMO_ASSERT( &m_mappers[m_row.at(r)] == &m_mappers[m_col.at(c)], "Error");
746
747 for (index_t i = 0; i != numActive; ++i)
748 {
749 const index_t ii = m_rstr.at(r) + actives.at(i); // N_i
750
751 for (index_t j = 0; j != numActive; ++j)
752 {
753 const index_t jj = m_cstr.at(c) + actives.at(j); // N_j
754
755 // If matrix is symmetric, we store only lower
756 // triangular part
757 if ( (!symm) || jj <= ii )
758 m_matrix.coeffRef(ii, jj) += localMat(i, j);
759
760
761 }
762 }
763 }
764
765
778 void pushToMatrixAllFree(const gsMatrix<T> & localMat,
779 const gsMatrix<index_t> & actives_i,
780 const gsMatrix<index_t> & actives_j,
781 const size_t r = 0, const size_t c = 0)
782 {
783 const index_t numActive_i = actives_i.rows();
784 const index_t numActive_j = actives_j.rows();
785
786 for (index_t i = 0; i != numActive_i; ++i)
787 {
788 const index_t ii = m_rstr.at(r) + actives_i.at(i); // N_i
789
790 for (index_t j = 0; j != numActive_j; ++j)
791 {
792 const index_t jj = m_cstr.at(c) + actives_j.at(j); // N_j
793
794 // If matrix is symmetric, we store only lower
795 // triangular part
796 if ( (!symm) || jj <= ii )
797 m_matrix.coeffRef(ii, jj) += localMat(i, j);
798 }
799 }
800 }
801
816 void pushToMatrix(const gsMatrix<T> & localMat,
817 const std::vector<gsMatrix<index_t> >& actives_vec,
818 const std::vector<gsMatrix<T> > & eliminatedDofs,
819 const gsVector<index_t> & r_vec,
820 const gsVector<index_t> & c_vec)
821 {
822 int rstrLocal = 0;
823 int cstrLocal = 0;
824
825 for (index_t r_ind = 0; r_ind != r_vec.size(); ++r_ind) // for row-blocks
826 {
827 size_t r = r_vec(r_ind);
828 const gsDofMapper & rowMap = m_mappers[m_row.at(r)];
829 const index_t numActive_i = actives_vec[r].rows();
830
831 for (index_t c_ind = 0; c_ind != c_vec.size(); ++c_ind) // for col-blocks
832 {
833 size_t c = c_vec(c_ind);
834 const gsDofMapper & colMap = m_mappers[m_col.at(c)];
835 const index_t numActive_j = actives_vec[c].rows();
836 const gsMatrix<T> & eliminatedDofs_j = eliminatedDofs[c];
837
838 for (index_t i = 0; i != numActive_i; ++i) // N_i
839 {
840 const int ii = m_rstr.at(r) + actives_vec[r].at(i); // row index global matrix
841 const int iiLocal = rstrLocal + i; // row index local matrix
842
843 if ( rowMap.is_free_index(actives_vec[r].at(i)) )
844 {
845
846 for (index_t j = 0; j != numActive_j; ++j) // N_j
847 {
848 const int jj = m_cstr.at(c) + actives_vec[c].at(j); // column index global matrix
849 const int jjLocal = cstrLocal + j; // column index local matrix
850
851 if ( colMap.is_free_index(actives_vec[c].at(j)) )
852 {
853 // If matrix is symmetric, we store only lower
854 // triangular part
855 if ( (!symm) || jj <= ii )
856 m_matrix.coeffRef(ii, jj) += localMat(iiLocal, jjLocal);
857 }
858 else // Fixed DoF
859 {
860 m_rhs.row(ii).noalias() -= localMat(iiLocal, jjLocal) * eliminatedDofs_j.row( colMap.global_to_bindex(actives_vec[c].at(j)));
861 }
862 }
863 }
864 }
865 cstrLocal += numActive_j;
866 }
867 cstrLocal = 0;
868 rstrLocal += numActive_i;
869 }
870
871 }
872
873
874
875public: /* Add local contributions to system right-hand side */
876
884 void pushToRhs(const gsMatrix<T> & localRhs,
885 const gsMatrix<index_t> & actives,
886 const size_t r = 0)
887 {
888 const gsDofMapper & mapper = m_mappers[m_row.at(r)];
889 const index_t numActive = actives.rows();
890
891 for (index_t i = 0; i != numActive; ++i)
892 {
893 const index_t ii = m_rstr.at(r) + actives.at(i);
894 if ( mapper.is_free_index(actives.at(i)) )
895 {
896 m_rhs.row(ii) += localRhs.row(i);
897 }
898 }
899 }
900
908 void pushToRhsAllFree(const gsMatrix<T> & localRhs,
909 const gsMatrix<index_t> & actives,
910 const size_t r = 0)
911 {
912 const index_t numActive = actives.rows();
913
914 for (index_t i = 0; i != numActive; ++i)
915 {
916 const index_t ii = m_rstr.at(r) + actives.at(i);
917 m_rhs.row(ii) += localRhs.row(i);
918 }
919 }
920
921
930 void pushToRhs(const gsMatrix<T> & localRhs,
931 const std::vector<gsMatrix<index_t> >& actives_vec,
932 const gsVector<index_t> & r_vec)
933 {
934 int rstrLocal = 0;
935
936 for (index_t r_ind = 0; r_ind != r_vec.size(); ++r_ind) // for row-blocks
937 {
938 index_t r = r_vec(r_ind);
939 const gsDofMapper & rowMap = m_mappers[m_row.at(r)];
940 const index_t numActive_i = actives_vec[r].rows();
941
942
943 for (index_t i = 0; i != numActive_i; ++i) // N_i
944 {
945 const int ii = m_rstr.at(r) + actives_vec[r].at(i); // row index global matrix
946 const int iiLocal = rstrLocal + i; // row index local matrix
947
948 if ( rowMap.is_free_index(actives_vec[r].at(i)) )
949 {
950 m_rhs.row(ii) += localRhs.row(iiLocal);
951 }
952 }
953 rstrLocal += numActive_i;
954 }
955
956 }
957
958public: /* Add local contributions to system matrix and right-hand side */
959
972 void push(const gsMatrix<T> & localMat,
973 const gsMatrix<T> & localRhs,
974 const gsMatrix<index_t> & actives,
975 const gsMatrix<T> & eliminatedDofs,
976 const size_t r = 0, const size_t c = 0)
977 {
978 const index_t numActive = actives.rows();
979 const gsDofMapper & rowMap = m_mappers[m_row.at(r)];
980
981 GISMO_ASSERT( &rowMap == &m_mappers[m_col.at(c)], "Error");
982 GISMO_ASSERT( m_matrix.cols() == m_rhs.rows(), "gsSparseSystem is not allocated");
983 //Assert eliminatedDofs.rows() == rowMap.boundarySize()
984
985 for (index_t i = 0; i != numActive; ++i)
986 {
987 const int ii = m_rstr.at(r) + actives(i);
988 if ( rowMap.is_free_index(actives.at(i)) )
989 {
990 m_rhs.row(ii) += localRhs.row(i);
991
992 for (index_t j = 0; j < numActive; ++j)
993 {
994 const int jj = m_cstr.at(c) + actives(j);
995 if ( rowMap.is_free_index(actives.at(j)) )
996 {
997 // If matrix is symmetric, we store only lower
998 // triangular part
999 if ( (!symm) || jj <= ii )
1000 m_matrix.coeffRef(ii, jj) += localMat(i, j);
1001 }
1002 else // if ( mapper.is_boundary_index(jj) ) // Fixed DoF?
1003 {
1004 m_rhs.row(ii).noalias() -= localMat(i, j) *
1005 eliminatedDofs.row( rowMap.global_to_bindex(actives.at(j)) );
1006 }
1007 }
1008 }
1009 }
1010 }
1011
1026 void push(const gsMatrix<T> & localMat,
1027 const gsMatrix<T> & localRhs,
1028 const gsMatrix<index_t> & actives_i,
1029 const gsMatrix<index_t> & actives_j,
1030 const gsMatrix<T> & eliminatedDofs_j,
1031 const size_t r = 0, const size_t c = 0)
1032 {
1033 const index_t numActive_i = actives_i.rows();
1034 const index_t numActive_j = actives_j.rows();
1035 const gsDofMapper & rowMap = m_mappers[m_row.at(r)];
1036 const gsDofMapper & colMap = m_mappers[m_col.at(c)];
1037
1038 GISMO_ASSERT( m_matrix.cols() == m_rhs.rows(), "gsSparseSystem is not allocated");
1039 //Assert eliminatedDofs.rows() == rowMap.boundarySize()
1040
1041 for (index_t i = 0; i != numActive_i; ++i)
1042 {
1043 const int ii = m_rstr.at(r) + actives_i.at(i);
1044 if ( rowMap.is_free_index(actives_i.at(i)) )
1045 {
1046 m_rhs.row(ii) += localRhs.row(i);
1047
1048 for (index_t j = 0; j < numActive_j; ++j)
1049 {
1050 const int jj = m_cstr.at(c) + actives_j.at(j);
1051 if ( colMap.is_free_index(actives_j.at(j)) )
1052 {
1053 // If matrix is symmetric, we store only lower
1054 // triangular part
1055 if ( (!symm) || jj <= ii )
1056 m_matrix.coeffRef(ii, jj) += localMat(i, j);
1057 }
1058 else // if ( mapper.is_boundary_index(jj) ) // Fixed DoF?
1059 {
1060 m_rhs.row(ii).noalias() -= localMat(i, j) *
1061 eliminatedDofs_j.row( colMap.global_to_bindex(actives_j.at(j)) );
1062 }
1063 }
1064 }
1065 }
1066 }
1067
1068
1082 void pushAllFree(const gsMatrix<T> & localMat,
1083 const gsMatrix<T> & localRhs,
1084 const gsMatrix<index_t> & actives,
1085 const size_t r = 0, const size_t c = 0)
1086 {
1087 const index_t numActive = actives.rows();
1088 //GISMO_ASSERT( &m_mappers[m_row.at(r)] == &m_mappers[m_col.at(c)], "Error");
1089
1090 for (index_t j=0; j!=numActive; ++j)
1091 {
1092 const unsigned jj = m_cstr.at(c) + actives(j);
1093 m_rhs.row(jj) += localRhs.row(j);
1094 for (index_t i=0; i!=numActive; ++i)
1095 {
1096 const unsigned ii = m_rstr.at(r) + actives(i);
1097 // If matrix is symmetric, we store only lower
1098 // triangular part
1099 if ( (!symm) || jj <= ii )
1100 m_matrix( ii, jj ) += localMat(i,j);
1101 }
1102 }
1103 }
1104
1105
1122 void push(const gsMatrix<T> & localMat,
1123 const gsMatrix<T> & localRhs,
1124 const std::vector<gsMatrix<index_t> >& actives_vec,
1125 const std::vector<gsMatrix<T> > & eliminatedDofs,
1126 const gsVector<index_t> & r_vec,
1127 const gsVector<index_t> & c_vec)
1128 {
1129 int rstrLocal = 0;
1130 int cstrLocal = 0;
1131
1132 for (index_t r_ind = 0; r_ind != r_vec.size(); ++r_ind) // for row-blocks
1133 {
1134 index_t r = r_vec(r_ind);
1135 const gsDofMapper & rowMap = m_mappers[m_row.at(r)];
1136 const index_t numActive_i = actives_vec[r].rows();
1137
1138 for (index_t c_ind = 0; c_ind != c_vec.size(); ++c_ind) // for col-blocks
1139 {
1140 index_t c = c_vec(c_ind);
1141 const gsDofMapper & colMap = m_mappers[m_col.at(c)];
1142 const index_t numActive_j = actives_vec[c].rows();
1143 const gsMatrix<T> & eliminatedDofs_j = eliminatedDofs[c];
1144
1145 for (index_t i = 0; i != numActive_i; ++i) // N_i
1146 {
1147 const int ii = m_rstr.at(r) + actives_vec[r].at(i); // row index global matrix
1148 const int iiLocal = rstrLocal + i; // row index local matrix
1149
1150 if ( rowMap.is_free_index(actives_vec[r].at(i)) )
1151 {
1152 // rhs should not be pushed for each col-block (but only once)
1153 if(c_ind == 0)
1154 m_rhs.row(ii) += localRhs.row(iiLocal);
1155
1156 for (index_t j = 0; j != numActive_j; ++j) // N_j
1157 {
1158 const int jj = m_cstr.at(c) + actives_vec[c].at(j); // column index global matrix
1159 const int jjLocal = cstrLocal + j; // column index local matrix
1160
1161 if ( colMap.is_free_index(actives_vec[c].at(j)) )
1162 {
1163 // If matrix is symmetric, we store only lower
1164 // triangular part
1165 if ( (!symm) || jj <= ii )
1166 m_matrix.coeffRef(ii, jj) += localMat(iiLocal, jjLocal);
1167 }
1168 else // Fixed DoF
1169 {
1170 m_rhs.row(ii).noalias() -= localMat(iiLocal, jjLocal) * eliminatedDofs_j.row( colMap.global_to_bindex(actives_vec[c].at(j)));
1171 }
1172 }
1173 }
1174 }
1175 cstrLocal += numActive_j;
1176 }
1177 cstrLocal = 0;
1178 rstrLocal += numActive_i;
1179 }
1180
1181 }
1182
1183
1184 // AM: incomplete, not working yet
1201 void push(const gsMatrix<T> & localMat,
1202 const gsMatrix<T> & localRhs,
1203 const std::vector<gsMatrix<index_t> >& actives,
1204 const std::vector<gsMatrix<T> > & eliminatedDofs)
1205 {
1206 GISMO_ASSERT( m_matrix.cols() == m_rhs.rows(), "gsSparseSystem is not allocated");
1207
1208 for (size_t r = 0; r != actives.size(); ++r) // for all row-blocks
1209 {
1210 const gsDofMapper & rowMap = m_mappers[m_row.at(r)];
1211 const index_t numRowActive = actives[r].rows();
1212
1213 for (size_t c = 0; c != actives.size(); ++c) // for all col-blocks
1214 {
1215 const gsMatrix<T> & fixedDofs = eliminatedDofs[m_col.at(c)];
1216
1217 for (index_t i = 0; i != numRowActive; ++i)
1218 {
1219 const int ii = m_rstr.at(r) + actives[r].at(i);
1220 if ( rowMap.is_free_index(actives[r].at(i)) )
1221 {
1222 m_rhs.row(ii) += localRhs.row(i + r * numRowActive); // + c *
1223 const index_t numColActive = actives[c].rows();
1224
1225 for (index_t j = 0; j < numColActive; ++j)
1226 {
1227 const int jj = m_cstr.at(c) + actives[c].at(j);
1228 if ( rowMap.is_free_index(actives[c].at(j)) )
1229 {
1230 // If matrix is symmetric, we store only lower
1231 // triangular part
1232 if ( (!symm) || jj <= ii )
1233 m_matrix.coeffRef(ii, jj) += localMat(i + r * numRowActive,
1234 j + c * numRowActive); // + c * ..
1235 }
1236 else // if ( mapper.is_boundary_index(jj) ) // Fixed DoF?
1237 {
1238 m_rhs.at(ii) -= localMat(i + r * numRowActive, j + c * numRowActive) * // + c *..
1239 fixedDofs.coeff( rowMap.global_to_bindex(actives[c].at(j)), 0 );
1240 }
1241 }
1242 }
1243 }
1244 }
1245 }
1246 }
1247
1264 void push(const std::vector<gsMatrix<T> > &,
1265 const std::vector<gsMatrix<T> > &,
1266 const std::vector<gsMatrix<index_t> > &,
1267 const std::vector<gsMatrix<T> > &,
1268 const gsVector<index_t> &, const gsVector<index_t> &)
1269 {
1270 GISMO_NO_IMPLEMENTATION
1271 }
1272
1273
1274 GISMO_DEPRECATED
1287 void pushToMatrix(const gsMatrix<T> & localMat,
1288 const gsMatrix<index_t> & actives,
1289 const size_t r = 0, const size_t c = 0)
1290 {
1291 const index_t numActive = actives.rows();
1292 const gsDofMapper & rowMap = m_mappers[m_row.at(r)];
1293 GISMO_ASSERT( &rowMap == &m_mappers[m_col.at(c)], "Error");
1294
1295 for (index_t i = 0; i != numActive; ++i)
1296 {
1297 const int ii = m_rstr.at(r) + actives.at(i); // N_i
1298
1299 if ( rowMap.is_free_index(actives.at(i)) )
1300 {
1301 for (index_t j = 0; j != numActive; ++j)
1302 {
1303 const int jj = m_cstr.at(c) + actives.at(j); // N_j
1304
1305 if ( rowMap.is_free_index( actives.at(j)) )
1306 // If matrix is symmetric, we store only lower
1307 // triangular part
1308 if ( (!symm) || jj <= ii )
1309 m_matrix.coeffRef(ii, jj) += localMat(i, j);
1310 }
1311 }
1312 }
1313 }
1314
1315 GISMO_DEPRECATED
1329 void pushToMatrix(const gsMatrix<T> & localMat,
1330 const gsMatrix<index_t> & actives_i,
1331 const gsMatrix<index_t> & actives_j,
1332 const size_t r = 0, const size_t c = 0)
1333 {
1334 const index_t numActive_i = actives_i.rows();
1335 const index_t numActive_j = actives_j.rows();
1336
1337 const gsDofMapper & rowMap = m_mappers[m_row.at(r)];
1338 const gsDofMapper & colMap = m_mappers[m_col.at(c)];
1339
1340 for (index_t i = 0; i != numActive_i; ++i)
1341 {
1342 const int ii = m_rstr.at(r) + actives_i.at(i); // N_i
1343
1344 if ( rowMap.is_free_index(actives_i.at(i)) )
1345 {
1346 for (index_t j = 0; j != numActive_j; ++j)
1347 {
1348 const int jj = m_cstr.at(c) + actives_j.at(j); // N_j
1349
1350 if ( colMap.is_free_index(actives_j.at(j)) )
1351 // If matrix is symmetric, we store only lower
1352 // triangular part
1353 if ( (!symm) || jj <= ii )
1354 m_matrix.coeffRef(ii, jj) += localMat(i, j);
1355 }
1356 }
1357 }
1358 }
1359
1360 GISMO_DEPRECATED
1373 void push(const gsMatrix<T> & localMat,
1374 const gsMatrix<T> & localRhs,
1375 const gsMatrix<index_t> & actives,
1376 const size_t r = 0, const size_t c = 0)
1377 {
1378 GISMO_ASSERT( m_matrix.cols() == m_rhs.rows(), "gsSparseSystem is not allocated");
1379
1380 const index_t numActive = actives.rows();
1381 const gsDofMapper & rowMap = m_mappers[m_row.at(r)];
1382 GISMO_ASSERT( &rowMap == &m_mappers[m_col.at(c)], "Error");
1383
1384 for (index_t i = 0; i != numActive; ++i)
1385 {
1386 const int ii = m_rstr.at(r) + actives(i);
1387 if ( rowMap.is_free_index(actives(i)) )
1388 {
1389 m_rhs.row(ii) += localRhs.row(i);
1390
1391 for (index_t j = 0; j != numActive; ++j)
1392 {
1393 const int jj = m_cstr.at(c) + actives(j);
1394 if ( rowMap.is_free_index(actives(j)) )
1395 {
1396 // If matrix is symmetric, we store only lower
1397 // triangular part
1398 if ( (!symm) || jj <= ii )
1399 m_matrix.coeffRef(ii, jj) += localMat(i, j);
1400 }
1401 }
1402 }
1403 }
1404 }
1405
1406
1407public:
1408 void pushSparse(const gsSparseMatrix<T> & localMat,
1409 const gsMatrix<T> & localRhs,
1410 const gsMatrix<index_t> & actives_i,
1411 const gsMatrix<index_t> & actives_j,
1412 const gsMatrix<T> & eliminatedDofs_j,
1413 const size_t r = 0, const size_t c = 0)
1414 {
1415 const gsDofMapper & rowMap = m_mappers[m_row.at(r)];
1416 const gsDofMapper & colMap = m_mappers[m_col.at(c)];
1417
1418 GISMO_ASSERT( m_matrix.cols() == m_rhs.rows(), "gsSparseSystem is not allocated");
1419 //Assert eliminatedDofs.rows() == rowMap.boundarySize()
1420
1421 for (index_t i = 0; i<localMat.outerSize(); ++i)
1422 for (typename gsSparseMatrix<T>::iterator it(localMat,i); it; ++it)
1423 {
1424 const int ii = m_rstr.at(r) + actives_i.at(it.row());
1425 if ( rowMap.is_free_index(actives_i.at(it.row())) )
1426 {
1427 const int jj = m_cstr.at(c) + actives_j.at(it.col());
1428 if ( colMap.is_free_index(actives_j.at(it.col())) )
1429 {
1430 // If matrix is symmetric, we store only lower
1431 // triangular part
1432 if ( (!symm) || jj <= ii )
1433 m_matrix.coeffRef(ii, jj) += it.value();
1434 }
1435 else // if ( mapper.is_boundary_index(jj) ) // Fixed DoF?
1436 {
1437 m_rhs.row(ii).noalias() -= it.value() *
1438 eliminatedDofs_j.row( colMap.global_to_bindex(jj) );
1439 }
1440 }
1441 }
1442
1443 for(index_t i=0; i<actives_i.rows();++i)
1444 {
1445 const int ii = m_rstr.at(r) + actives_i.at(i);
1446 if ( rowMap.is_free_index(actives_i.at(i)) )
1447 m_rhs.row(ii) += localRhs.row(i);
1448 }
1449 }
1450
1451}; // class gsSparseSystem
1452
1453
1454template<class T>
1455std::ostream &operator<<(std::ostream &os, gsSparseSystem<T> & ss)
1456{
1457 os << ss.blockView();
1458 os << "Right-hand side: "
1459 << ss.rhs().rows() <<"x"<< ss.rhs().cols() <<"\n";
1460 return os;
1461}
1462
1463
1464} // namespace gismo
gismo::gsBasis
A basis represents a family of scalar basis functions defined over a common parameter domain.
Definition gsBasis.h:79
gismo::gsBasis::degree
virtual short_t degree(short_t i) const
Degree with respect to the i-th variable. If the basis is a tensor product of (piecewise) polynomial ...
Definition gsBasis.hpp:699
gismo::gsDofMapper
Maintains a mapping from patch-local dofs to global dof indices and allows the elimination of individ...
Definition gsDofMapper.h:69
gismo::gsDofMapper::is_free_index
bool is_free_index(index_t gl) const
Returns true if global dof gl is not eliminated.
Definition gsDofMapper.h:376
gismo::gsDofMapper::global_to_bindex
index_t global_to_bindex(index_t gl) const
Returns the boundary index of global dof gl.
Definition gsDofMapper.h:364
gismo::gsMatrixBlockView
Represents a block-view of the given matrix.
Definition gsMatrixBlockView.h:32
gismo::gsMatrix
A matrix with arbitrary coefficient type and fixed or dynamic size.
Definition gsMatrix.h:41
gismo::gsMatrix::at
T at(index_t i) const
Returns the i-th element of the vectorization of the matrix.
Definition gsMatrix.h:211
gismo::gsMatrix::blockView
BlockView blockView(const gsVector< index_t > &rowSizes, const gsVector< index_t > &colSizes)
Return a block view of the matrix with rowSizes and colSizes.
Definition gsMatrix.h:381
gismo::gsMultiBasis
Holds a set of patch-wise bases and their topology information.
Definition gsMultiBasis.h:37
gismo::gsOptionList
Class which holds a list of parameters/options, and provides easy access to them.
Definition gsOptionList.h:33
gismo::gsOptionList::getInt
const index_t & getInt(const std::string &label) const
Reads value for option label from options.
Definition gsOptionList.cpp:37
gismo::gsOptionList::getReal
Real getReal(const std::string &label) const
Reads value for option label from options.
Definition gsOptionList.cpp:44
gismo::gsSparseMatrixIter
Iterator over the non-zero entries of a sparse matrix.
Definition gsSparseMatrix.h:74
gismo::gsSparseMatrix
Sparse matrix class, based on gsEigen::SparseMatrix.
Definition gsSparseMatrix.h:139
gismo::gsSparseMatrix::blockView
BlockView blockView(const gsVector< index_t > &rowSizes, const gsVector< index_t > &colSizes)
Return a block view of the matrix with rowSizes and colSizes.
Definition gsSparseMatrix.h:336
gismo::gsSparseSystem
A sparse linear system indexed by sets of degrees of freedom.
Definition gsSparseSystem.h:30
gismo::gsSparseSystem::m_rhs
gsMatrix< T > m_rhs
the right hand side of the system
Definition gsSparseSystem.h:46
gismo::gsSparseSystem::colMapper
gsDofMapper & colMapper(const index_t c)
colMapper returns the dofMapper for column block c
Definition gsSparseSystem.h:508
gismo::gsSparseSystem::pushToMatrixAllFree
void pushToMatrixAllFree(const gsMatrix< T > &localMat, const gsMatrix< index_t > &actives_i, const gsMatrix< index_t > &actives_j, const size_t r=0, const size_t c=0)
pushToMatrixAllFree pushes the local matrix for an element to the global system,
Definition gsSparseSystem.h:778
gismo::gsSparseSystem::rows
index_t rows() const
the number of matrix rows
Definition gsSparseSystem.h:383
gismo::gsSparseSystem::numUnknowns
index_t numUnknowns() const
returns the number of unknowns
Definition gsSparseSystem.h:524
gismo::gsSparseSystem::rhsRows
index_t rhsRows() const
the rows of the right-hand side vector
Definition gsSparseSystem.h:386
gismo::gsSparseSystem::mapRowIndices
void mapRowIndices(const gsMatrix< index_t > &actives, const index_t patchIndex, gsMatrix< index_t > &result, const index_t r=0) const
mapRowIndices Maps a set of basis indices by the corresponding dofMapper.
Definition gsSparseSystem.h:567
gismo::gsSparseSystem::rhs
const gsMatrix< T > & rhs() const
Access the right hand side.
Definition gsSparseSystem.h:402
gismo::gsSparseSystem::gsSparseSystem
gsSparseSystem(DofMappers &mappers, const index_t rows, const index_t cols)
gsSparseSystem Constructor for the sparse system, with given number of row and column blocks....
Definition gsSparseSystem.h:205
gismo::gsSparseSystem::pushToMatrix
void pushToMatrix(const gsMatrix< T > &localMat, const gsMatrix< index_t > &actives, const gsMatrix< T > &eliminatedDofs, const size_t r=0, const size_t c=0)
pushToMatrix pushes the local matrix for an element to the global system,
Definition gsSparseSystem.h:638
gismo::gsSparseSystem::pushAllFree
void pushAllFree(const gsMatrix< T > &localMat, const gsMatrix< T > &localRhs, const gsMatrix< index_t > &actives, const size_t r=0, const size_t c=0)
pushAllFree pushes the local system matrix and rhs for an element to the global system,
Definition gsSparseSystem.h:1082
gismo::gsSparseSystem::reserve
void reserve(const gsMultiBasis< T > &mb, const gsOptionList &opt, const index_t numRhs)
Reserves the memory for the sparse matrix and the rhs, based on the polynomial degree of the first ba...
Definition gsSparseSystem.h:350
gismo::gsSparseSystem::push
void push(const gsMatrix< T > &localMat, const gsMatrix< T > &localRhs, const gsMatrix< index_t > &actives_i, const gsMatrix< index_t > &actives_j, const gsMatrix< T > &eliminatedDofs_j, const size_t r=0, const size_t c=0)
push pushes the local system matrix and rhs for an element to the global system,
Definition gsSparseSystem.h:1026
gismo::gsSparseSystem::pushToRhsAllFree
void pushToRhsAllFree(const gsMatrix< T > &localRhs, const gsMatrix< index_t > &actives, const size_t r=0)
pushToRhsAllFree pushes the local rhs for an element to the global system
Definition gsSparseSystem.h:908
gismo::gsSparseSystem::dofMappers
const DofMappers & dofMappers() const
returns all dof Mappers.
Definition gsSparseSystem.h:528
gismo::gsSparseSystem::colMapper
const gsDofMapper & colMapper(const index_t c) const
colMapper returns the dofMapper for column block c
Definition gsSparseSystem.h:498
gismo::gsSparseSystem::setZero
void setZero()
set everything to zero
Definition gsSparseSystem.h:373
gismo::gsSparseSystem::pushToMatrix
void pushToMatrix(const gsMatrix< T > &localMat, const std::vector< gsMatrix< index_t > > &actives_vec, const std::vector< gsMatrix< T > > &eliminatedDofs, const gsVector< index_t > &r_vec, const gsVector< index_t > &c_vec)
pushToMatrix pushes one local matrix consisting of several blocks corresponding to blocks of the glob...
Definition gsSparseSystem.h:816
gismo::gsSparseSystem::pushToMatrix
GISMO_DEPRECATED void pushToMatrix(const gsMatrix< T > &localMat, const gsMatrix< index_t > &actives_i, const gsMatrix< index_t > &actives_j, const size_t r=0, const size_t c=0)
pushToMatrix pushes the local matrix for an element to the global system,
Definition gsSparseSystem.h:1329
gismo::gsSparseSystem::mapToGlobalColIndex
void mapToGlobalColIndex(const index_t active, const index_t patchIndex, index_t &result, const index_t c=0) const
mapToGlobalColIndex Maps a single basis index to the final position in the system,...
Definition gsSparseSystem.h:616
gismo::gsSparseSystem::mapToGlobalRowIndex
void mapToGlobalRowIndex(const index_t active, const index_t patchIndex, index_t &result, const index_t r=0) const
mapToGlobalRowIndex Maps a single basis index to the final position in the system,...
Definition gsSparseSystem.h:600
gismo::gsSparseSystem::blockView
matBlockView blockView()
returns a block view of the matrix, easy way to extract single blocks
Definition gsSparseSystem.h:410
gismo::gsSparseSystem::pushToRhs
void pushToRhs(const gsMatrix< T > &localRhs, const gsMatrix< index_t > &actives, const size_t r=0)
pushToRhs pushes the local rhs for an element to the global system
Definition gsSparseSystem.h:884
gismo::gsSparseSystem::m_rstr
gsVector< index_t > m_rstr
strides for the row blocks (shifting of mapped indices). The mapper do not have this information anym...
Definition gsSparseSystem.h:68
gismo::gsSparseSystem::symmetry
bool symmetry() const
returns true if only half of the matrix is stored, due to its symmetry
Definition gsSparseSystem.h:538
gismo::gsSparseSystem::gsSparseSystem
gsSparseSystem(DofMappers &mappers, const gsVector< index_t > &rowInd, const gsVector< index_t > &colInd, const gsVector< index_t > &colvar)
gsSparseSystem The constructor for the sparse system given the set of mappers, may be different for e...
Definition gsSparseSystem.h:275
gismo::gsSparseSystem::pushToMatrix
void pushToMatrix(const gsMatrix< T > &localMat, const gsMatrix< index_t > &actives_i, const gsMatrix< index_t > &actives_j, const gsMatrix< T > &eliminatedDofs_j, const size_t r=0, const size_t c=0)
pushToMatrix pushes the local matrix for an element to the global system,
Definition gsSparseSystem.h:689
gismo::gsSparseSystem::m_mappers
DofMappers m_mappers
the mappers for the different blocks. The assignment of the mappers to the column (and row) blocks is...
Definition gsSparseSystem.h:56
gismo::gsSparseSystem::m_cvar
gsVector< index_t > m_cvar
map between column blocks and used multibasis (which unknown/component correlate to which multibasis)...
Definition gsSparseSystem.h:77
gismo::gsSparseSystem::gsSparseSystem
gsSparseSystem(DofMappers &mappers, const gsVector< index_t > &dims)
gsSparseSystem Constructor of a sparse System specified by the number of unknows for each block,...
Definition gsSparseSystem.h:139
gismo::gsSparseSystem::gsSparseSystem
gsSparseSystem(gsDofMapper &mapper)
gsSparseSystem Constuctor for the sparse System only considering one block (scalar equation),...
Definition gsSparseSystem.h:91
gismo::gsSparseSystem::numColNz
index_t numColNz(const gsMultiBasis< T > &mb, const gsOptionList &opt) const
Provides an estimation of the number of non-zero matrix entries per column. This value can be used fo...
Definition gsSparseSystem.h:359
gismo::gsSparseSystem::rhsCols
index_t rhsCols() const
the number of right-hand side vectors
Definition gsSparseSystem.h:389
gismo::gsSparseSystem::m_cstr
gsVector< index_t > m_cstr
strides for the column blocks (shifting of mapped indices). The mapper do not have this information a...
Definition gsSparseSystem.h:72
gismo::gsSparseSystem::pushToMatrix
GISMO_DEPRECATED void pushToMatrix(const gsMatrix< T > &localMat, const gsMatrix< index_t > &actives, const size_t r=0, const size_t c=0)
pushToMatrix pushes the local matrix for an element to the global system,
Definition gsSparseSystem.h:1287
gismo::gsSparseSystem::pushToMatrixAllFree
void pushToMatrixAllFree(const gsMatrix< T > &localMat, const gsMatrix< index_t > &actives, const size_t r=0, const size_t c=0)
pushToMatrixAllFree pushes the local matrix for an element to the global system,
Definition gsSparseSystem.h:740
gismo::gsSparseSystem::reserve
void reserve(const index_t nz, const index_t numRhs)
reserve reserves the memory for the sparse matrix and the rhs.
Definition gsSparseSystem.h:327
gismo::gsSparseSystem::rhs
gsMatrix< T > & rhs()
Access the right hand side.
Definition gsSparseSystem.h:406
gismo::gsSparseSystem::matrix
gsSparseMatrix< T > & matrix()
Access the system Matrix.
Definition gsSparseSystem.h:398
gismo::gsSparseSystem::m_row
gsVector< index_t > m_row
map between row blocks and index of m_mappers, i.e. row block i is described by m_mappers[m_row[i]].
Definition gsSparseSystem.h:60
gismo::gsSparseSystem::push
void push(const gsMatrix< T > &localMat, const gsMatrix< T > &localRhs, const gsMatrix< index_t > &actives, const gsMatrix< T > &eliminatedDofs, const size_t r=0, const size_t c=0)
push pushes the local system matrix and rhs for an element to the global system,
Definition gsSparseSystem.h:972
gismo::gsSparseSystem::push
void push(const gsMatrix< T > &localMat, const gsMatrix< T > &localRhs, const std::vector< gsMatrix< index_t > > &actives_vec, const std::vector< gsMatrix< T > > &eliminatedDofs, const gsVector< index_t > &r_vec, const gsVector< index_t > &c_vec)
pushToMatrix pushes one local matrix and rhs consisting of several blocks corresponding to blocks of ...
Definition gsSparseSystem.h:1122
gismo::gsSparseSystem::rowMapper
const gsDofMapper & rowMapper(const index_t r) const
rowMapper returns the dofMapper for row block r
Definition gsSparseSystem.h:482
gismo::gsSparseSystem::push
void push(const gsMatrix< T > &localMat, const gsMatrix< T > &localRhs, const std::vector< gsMatrix< index_t > > &actives, const std::vector< gsMatrix< T > > &eliminatedDofs)
push pushes one local system matrix and one local rhs for an element to the global system for several...
Definition gsSparseSystem.h:1201
gismo::gsSparseSystem::blockView
matBlockView blockView(size_t numRowBlocksNew, size_t numColBlocksNew, const gsVector< index_t > &rowBlocksNew, const gsVector< index_t > &colBlocksNew)
returns a block view of the matrix, where you can choose which blocks (of gsSparseSystem) should be c...
Definition gsSparseSystem.h:438
gismo::gsSparseSystem::swap
void swap(gsSparseSystem &other)
swap swaps the content of the Sparse System with the other given one
Definition gsSparseSystem.h:309
gismo::gsSparseSystem::blockViewRhs
rhsBlockView blockViewRhs(size_t numRowBlocksNew, const gsVector< index_t > &rowBlocksNew)
returns a block view of the rhs, where you can choose which blocks (of gsSparseSystem) should be comb...
Definition gsSparseSystem.h:459
gismo::gsSparseSystem::mapColIndices
void mapColIndices(const gsMatrix< index_t > &actives, const index_t patchIndex, gsMatrix< index_t > &result, const index_t c=0) const
mapColIndices Maps a set of basis indices by the corresponding dofMapper.
Definition gsSparseSystem.h:584
gismo::gsSparseSystem::pushToRhs
void pushToRhs(const gsMatrix< T > &localRhs, const std::vector< gsMatrix< index_t > > &actives_vec, const gsVector< index_t > &r_vec)
pushToRhs pushes one local rhs consisting of several blocks corresponding to blocks of the global sys...
Definition gsSparseSystem.h:930
gismo::gsSparseSystem::cols
index_t cols() const
the number of matrix columns
Definition gsSparseSystem.h:380
gismo::gsSparseSystem::m_col
gsVector< index_t > m_col
map between column blocks and index of m_mappers i.e. col block j is described by m_mappers[m_col[j]]...
Definition gsSparseSystem.h:64
gismo::gsSparseSystem::rowMapper
gsDofMapper & rowMapper(const index_t r)
rowMapper returns the dofMapper for row block r
Definition gsSparseSystem.h:490
gismo::gsSparseSystem::unkSize
index_t unkSize(const index_t unk) const
return the number of components for the given component
Definition gsSparseSystem.h:520
gismo::gsSparseSystem::m_matrix
gsSparseMatrix< T > m_matrix
the system matrix
Definition gsSparseSystem.h:43
gismo::gsSparseSystem::push
GISMO_DEPRECATED void push(const gsMatrix< T > &localMat, const gsMatrix< T > &localRhs, const gsMatrix< index_t > &actives, const size_t r=0, const size_t c=0)
push pushes the local system matrix and rhs for an element to the global system,
Definition gsSparseSystem.h:1373
gismo::gsSparseSystem::numColBlocks
index_t numColBlocks() const
returns the number of column blocks
Definition gsSparseSystem.h:472
gismo::gsSparseSystem::colBasis
index_t colBasis(const index_t c) const
colBasis returns the index of the Basis used for column block c
Definition gsSparseSystem.h:516
gismo::gsSparseSystem::matrix
const gsSparseMatrix< T > & matrix() const
Access the system Matrix.
Definition gsSparseSystem.h:394
gismo::gsSparseSystem::push
void push(const std::vector< gsMatrix< T > > &, const std::vector< gsMatrix< T > > &, const std::vector< gsMatrix< index_t > > &, const std::vector< gsMatrix< T > > &, const gsVector< index_t > &, const gsVector< index_t > &)
push pushes several local system matrices and local rhs for an element to the global system for sever...
Definition gsSparseSystem.h:1264
gismo::gsSparseSystem::numRowBlocks
index_t numRowBlocks() const
returns the number of row blocks
Definition gsSparseSystem.h:475
gismo::gsSparseSystem::initialized
bool initialized() const
checks if the system was initialized
Definition gsSparseSystem.h:532
gismo::gsStdVectorRef
Simple wrapper class for a vector of objects.
Definition gsStdVectorRef.h:28
gismo::gsVector
A vector with arbitrary coefficient type and fixed or dynamic size.
Definition gsVector.h:37
gismo::gsVector::at
T at(index_t i) const
Returns the i-th element of the vector.
Definition gsVector.h:177
short_t
#define short_t
Definition gsConfig.h:35
index_t
#define index_t
Definition gsConfig.h:32
GISMO_NO_IMPLEMENTATION
#define GISMO_NO_IMPLEMENTATION
Definition gsDebug.h:129
GISMO_ERROR
#define GISMO_ERROR(message)
Definition gsDebug.h:118
GISMO_ENSURE
#define GISMO_ENSURE(cond, message)
Definition gsDebug.h:102
GISMO_ASSERT
#define GISMO_ASSERT(cond, message)
Definition gsDebug.h:89
gsStdVectorRef.h
Small wrapper for std::vector.
gismo::expr::operator<<
std::ostream & operator<<(std::ostream &os, const _expr< E > &b)
Stream operator for expressions.
Definition gsExpressions.h:382
gismo
The G+Smo namespace, containing all definitions for the library.