G+Smo  24.08.0
Geometry + Simulation Modules
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
gsSparseSystem.h
Go to the documentation of this file.
1 
15 #pragma once
16 
17 #include <gsCore/gsStdVectorRef.h>
18 
19 namespace gismo
20 {
21 
28 template<typename T, bool symm>
29 class gsSparseSystem
30 {
31 protected:
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 
40 protected:
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 
81 public:
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 
391 public: /* 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 
556 public: /* 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 
624 public: /* 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 
875 public: /* 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 
958 public: /* 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 
1407 public:
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 
1454 template<class T>
1455 std::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::gsSparseSystem::blockView
matBlockView blockView()
returns a block view of the matrix, easy way to extract single blocks
Definition: gsSparseSystem.h:410
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::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::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::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::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::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:298
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_NO_IMPLEMENTATION
#define GISMO_NO_IMPLEMENTATION
Definition: gsDebug.h:129
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::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::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::gsStdVectorRef
Simple wrapper class for a vector of objects.
Definition: gsStdVectorRef.h:27
gismo::gsSparseSystem::rhsCols
index_t rhsCols() const
the number of right-hand side vectors
Definition: gsSparseSystem.h:389
short_t
#define short_t
Definition: gsConfig.h:35
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::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(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::cols
index_t cols() const
the number of matrix columns
Definition: gsSparseSystem.h:380
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::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, i.e. including the right shift due to the block structure
Definition: gsSparseSystem.h:616
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:650
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::m_matrix
gsSparseMatrix< T > m_matrix
the system matrix
Definition: gsSparseSystem.h:43
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::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::numColBlocks
index_t numColBlocks() const
returns the number of column blocks
Definition: gsSparseSystem.h:472
index_t
#define index_t
Definition: gsConfig.h:32
GISMO_ENSURE
#define GISMO_ENSURE(cond, message)
Definition: gsDebug.h:102
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::gsOptionList::getInt
const index_t & getInt(const std::string &label) const
Reads value for option label from options.
Definition: gsOptionList.cpp:37
gismo::gsMatrix
A matrix with arbitrary coefficient type and fixed or dynamic size.
Definition: gsMatrix.h:38
gismo::gsSparseSystem::rows
index_t rows() const
the number of matrix rows
Definition: gsSparseSystem.h:383
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_ASSERT
#define GISMO_ASSERT(cond, message)
Definition: gsDebug.h:89
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::gsDofMapper
Maintains a mapping from patch-local dofs to global dof indices and allows the elimination of individ...
Definition: gsDofMapper.h:68
gismo::gsSparseMatrix
Sparse matrix class, based on gsEigen::SparseMatrix.
Definition: gsSparseMatrix.h:140
gismo::gsSparseSystem::gsSparseSystem
gsSparseSystem(gsDofMapper &mapper)
gsSparseSystem Constuctor for the sparse System only considering one block (scalar equation)...
Definition: gsSparseSystem.h:91
gismo::gsSparseSystem::dofMappers
const DofMappers & dofMappers() const
returns all dof Mappers.
Definition: gsSparseSystem.h:528
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::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::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::gsSparseSystem::rowMapper
gsDofMapper & rowMapper(const index_t r)
rowMapper returns the dofMapper for row block r
Definition: gsSparseSystem.h:490
gismo::gsSparseSystem::initialized
bool initialized() const
checks if the system was initialized
Definition: gsSparseSystem.h:532
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::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::matrix
const gsSparseMatrix< T > & matrix() const
Access the system Matrix.
Definition: gsSparseSystem.h:394
gismo::gsMultiBasis
Holds a set of patch-wise bases and their topology information.
Definition: gsMultiBasis.h:36
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
gsStdVectorRef.h
Small wrapper for std::vector.
gismo::gsSparseSystem::colMapper
gsDofMapper & colMapper(const index_t c)
colMapper returns the dofMapper for column block c
Definition: gsSparseSystem.h:508
gismo::gsVector::at
T at(index_t i) const
Returns the i-th element of the vector.
Definition: gsVector.h:177
gismo::gsSparseSystem::m_rhs
gsMatrix< T > m_rhs
the right hand side of the system
Definition: gsSparseSystem.h:46
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::rowMapper
const gsDofMapper & rowMapper(const index_t r) const
rowMapper returns the dofMapper for row block r
Definition: gsSparseSystem.h:482
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::matrix
gsSparseMatrix< T > & matrix()
Access the system Matrix.
Definition: gsSparseSystem.h:398
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::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::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::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::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::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::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_ERROR
#define GISMO_ERROR(message)
Definition: gsDebug.h:118
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::numUnknowns
index_t numUnknowns() const
returns the number of unknowns
Definition: gsSparseSystem.h:524
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::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, i.e. including the right shift due to the block structure
Definition: gsSparseSystem.h:600
gismo::gsOptionList
Class which holds a list of parameters/options, and provides easy access to them. ...
Definition: gsOptionList.h:32
gismo::gsSparseSystem::rhs
gsMatrix< T > & rhs()
Access the right hand side.
Definition: gsSparseSystem.h:406
gismo::gsSparseSystem::rhs
const gsMatrix< T > & rhs() const
Access the right hand side.
Definition: gsSparseSystem.h:402
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::swap
void swap(gsSparseSystem &other)
swap swaps the content of the Sparse System with the other given one
Definition: gsSparseSystem.h:309
gismo::gsSparseSystem
A sparse linear system indexed by sets of degrees of freedom.
Definition: gsSparseSystem.h:29
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::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::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::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::gsBasis
A basis represents a family of scalar basis functions defined over a common parameter domain...
Definition: gsBasis.h:78
gismo::gsOptionList::getReal
Real getReal(const std::string &label) const
Reads value for option label from options.
Definition: gsOptionList.cpp:44
gismo::gsSparseSystem::rhsRows
index_t rhsRows() const
the rows of the right-hand side vector
Definition: gsSparseSystem.h:386
gismo::gsSparseSystem::numRowBlocks
index_t numRowBlocks() const
returns the number of row blocks
Definition: gsSparseSystem.h:475