OpenFCST: The open-source Fuel Cell Simulation Toolbox
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
mesh_worker_assembler.h
Go to the documentation of this file.
1 //---------------------------------------------------------------------------
2 // $Id: mesh_worker_assembler.h 1348 2013-08-16 00:45:07Z secanell $
3 //
4 // Copyright (C) 2006, 2007, 2008, 2009 by Guido Kanschat
5 //
6 // This file is subject to QPL and may not be distributed
7 // without copyright and license information. Please refer
8 // to the file deal.II/doc/license.html for the text and
9 // further information on this license.
10 //
11 //---------------------------------------------------------------------------
12 
13 #ifndef __deal2__mesh_worker_assembler_h
14 #define __deal2__mesh_worker_assembler_h
15 
16 namespace MeshWorker
17 {
18 
69  namespace Assembler
70  {
79  template <typename number = double>
80  class Functional
81  {
82  public:
83  void initialize(unsigned int n);
89  template<int dim>
90  void assemble(const InfoObjects::DoFInfo<dim>& info);
91 
97  template<int dim>
98  void assemble(const InfoObjects::DoFInfo<dim>& info1,
99  const InfoObjects::DoFInfo<dim>& info2);
100 
104  number operator() (unsigned int i) const;
105 
106  std::vector<number> J1;
107  private:
112  std::vector<double> results;
113  };
114 
121  template <typename number = double>
122  class CellsAndFaces
123  {
124  public:
125  void initialize(FEVectors& residuals, bool separate_faces = true);
131  template<int dim>
132  void assemble(const InfoObjects::DoFInfo<dim>& info);
133 
139  template<int dim>
140  void assemble(const InfoObjects::DoFInfo<dim>& info1,
141  const InfoObjects::DoFInfo<dim>& info2);
142 
146  number operator() (unsigned int i) const;
147 
148  std::vector<number> J1;
149  std::vector<number> J2;
150  private:
151  FEVectors results;
152  bool separate_faces;
153  };
154 
155 
176  template <typename number = double>
177  class ResidualLocalBlocksToGlobalBlocks
178  {
179  public:
195  void initialize(const BlockInfo& block_info,
196  FEVectors& residuals);
197 
203  template<int dim>
204  void assemble(const InfoObjects::DoFInfo<dim>& info);
205 
211  template<int dim>
212  void assemble(const InfoObjects::DoFInfo<dim>& info1,
213  const InfoObjects::DoFInfo<dim>& info2);
214 
220  std::vector<BlockVector<number> > R1;
221 
226  std::vector<BlockVector<number> > R2;
227  private:
232  void assemble( FEVector& global,
233  const BlockVector<number>& local,
234  const std::vector<unsigned int>& dof);
235 
240  BlockInfo block_info;
246  FEVectors residuals;
247  };
248 
255  template <typename number>
256  class LocalMatrixBlocks
257  {
263  void initialize_local(
264  MatrixBlock<FullMatrix<number> >& M,
265  unsigned int row, unsigned int col,
266  unsigned int n_rows, unsigned int n_cols);
267  public:
274  template <class MatrixPtr>
275  void initialize(const BlockInfo& block_info,
276  std::vector<MatrixPtr>& matrices);
277 
285  std::vector<MatrixBlock<FullMatrix<number> > > M11;
286 
293  std::vector<MatrixBlock<FullMatrix<number> > > M22;
294 
304  std::vector<MatrixBlock<FullMatrix<number> > > M12;
305 
315  std::vector<MatrixBlock<FullMatrix<number> > > M21;
316 
321  BlockInfo block_info;
322  };
323 
324 
325 
352  template <class MATRIX, typename number = double>
353  class MatrixLocalBlocksToGlobalBlocks
354  : public LocalMatrixBlocks<number>
355 
356  {
357  public:
359  typedef boost::shared_ptr<MatrixBlock<MATRIX> > MatrixPtr;
360 
368  MatrixLocalBlocksToGlobalBlocks(double threshold = 1.e-12);
369 
376  void initialize(const BlockInfo& block_info,
377  std::vector<MatrixPtr>& matrices);
378 
384  template<int dim>
385  void assemble(const InfoObjects::DoFInfo<dim>& info);
386 
392  template<int dim>
393  void assemble(const InfoObjects::DoFInfo<dim>& info1,
394  const InfoObjects::DoFInfo<dim>& info2);
395 
396  private:
401  void assemble(
402  MATRIX& global,
403  const FullMatrix<number>& local,
404  unsigned int block_row,
405  unsigned int block_col,
406  const std::vector<unsigned int>& dof1,
407  const std::vector<unsigned int>& dof2);
408 
414  std::vector<MatrixPtr> matrices;
415 
425  const double threshold;
426 
427  };
428 
452  template <class MATRIX, typename number = double>
453  class MGMatrixLocalBlocksToGlobalBlocks
454  : public LocalMatrixBlocks<number>
455  {
456  public:
458  typedef boost::shared_ptr<MatrixBlock<MGLevelObject<MATRIX> > > MatrixPtr;
459 
467  MGMatrixLocalBlocksToGlobalBlocks(double threshold = 1.e-12);
468 
475  void initialize(const BlockInfo& block_info,
476  std::vector<MatrixPtr>& matrices);
477 
489  void initialize_edge_flux(std::vector<MatrixPtr>& up, std::vector<MatrixPtr>& down);
490 
496  template<int dim>
497  void assemble(const InfoObjects::DoFInfo<dim>& info);
498 
504  template<int dim>
505  void assemble(const InfoObjects::DoFInfo<dim>& info1,
506  const InfoObjects::DoFInfo<dim>& info2);
507 
508  private:
513  void assemble(
514  MATRIX& global,
515  const FullMatrix<number>& local,
516  unsigned int block_row,
517  unsigned int block_col,
518  const std::vector<unsigned int>& dof1,
519  const std::vector<unsigned int>& dof2,
520  unsigned int level1,
521  unsigned int level2,
522  bool transpose = false);
523 
529  std::vector<MatrixPtr> matrices;
530 
536  std::vector<MatrixPtr> flux_down;
537 
543  std::vector<MatrixPtr> flux_up;
544 
554  const double threshold;
555 
556  };
557 
558 //----------------------------------------------------------------------//
559 
560  template <typename number>
561  inline void
562  Functional<number>::initialize(unsigned int n)
563  {
564  J1.resize(n);
565  results.resize(n);
566  }
567 
568 
569  template <typename number>
570  template<int dim>
571  inline void
572  Functional<number>::assemble(const InfoObjects::DoFInfo<dim>& info)
573  {
574  for (unsigned int i=0;i<results.size();++i)
575  {
576  results[i] += J1[i];
577  J1[i] = 0.;
578  }
579  }
580 
581 
582  template <typename number>
583  template<int dim>
584  inline void
585  Functional<number>::assemble(const InfoObjects::DoFInfo<dim>& info1,
586  const InfoObjects::DoFInfo<dim>& info2)
587  {
588  for (unsigned int i=0;i<results.size();++i)
589  {
590  results[i] += J1[i];
591  J1[i] = 0.;
592  }
593  }
594 
595 
596  template <typename number>
597  inline number
598  Functional<number>::operator() (unsigned int i) const
599  {
600  AssertIndexRange(i, results.size());
601  return results[i];
602  }
603 
604 //----------------------------------------------------------------------//
605 
606  template <typename number>
607  inline void
608  CellsAndFaces<number>::initialize(FEVectors& r, bool sep)
609  {
610  Assert(r.vector_name(0) == "cells", FEVectors::ExcNameMismatch(0, "cells"));
611  Assert(r.vector_name(1) == "faces", FEVectors::ExcNameMismatch(1, "faces"));
612  AssertDimension(r.vector(0).n_blocks(), r.vector(1).n_blocks());
613 
614  results = r;
615  separate_faces = sep;
616 
617  J1.resize(results.vector(0).n_blocks());
618  J2.resize(results.vector(0).n_blocks());
619  }
620 
621 
622  template <typename number>
623  template<int dim>
624  inline void
625  CellsAndFaces<number>::assemble(const InfoObjects::DoFInfo<dim>& info)
626  {
627  for (unsigned int i=0;i<J1.size();++i)
628  {
629  if (separate_faces &&
630  info.face_number != deal_II_numbers::invalid_unsigned_int)
631  results.vector(1).block(i)(info.face->user_index()) += J1[i];
632 
633  else
634  results.vector(0).block(i)(info.cell->user_index()) += J1[i];
635  J1[i] = 0.;
636  }
637  }
638 
639 
640  template <typename number>
641  template<int dim>
642  inline void
643  CellsAndFaces<number>::assemble(const InfoObjects::DoFInfo<dim>& info1,
644  const InfoObjects::DoFInfo<dim>& info2)
645  {
646  for (unsigned int i=0;i<J1.size();++i)
647  {
648  if (separate_faces)
649  {
650  results.vector(1).block(i)(info1.face->user_index()) += J1[i];
651  if (info2.face != info1.face)
652  results.vector(1).block(i)(info2.face->user_index()) += J1[i];
653  J1[i] = 0.;
654  }
655  else
656  {
657  results.vector(0).block(i)(info1.cell->user_index()) += J1[i];
658  results.vector(0).block(i)(info2.cell->user_index()) += J2[i];
659  J1[i] = 0.;
660  J2[i] = 0.;
661  }
662  }
663  }
664 
665 
666 //----------------------------------------------------------------------//
667 
668  template <typename number>
669  inline void
670  ResidualLocalBlocksToGlobalBlocks<number>::initialize(
671  const BlockInfo& bi,
672  FEVectors& m)
673  {
674  block_info = bi;
675  residuals = m;
676 
677  Assert(block_info.local.size() != 0, ExcNotInitialized());
678 
679  R1.resize(residuals.n_vectors());
680  R2.resize(residuals.n_vectors());
681  for (unsigned int i=0;i<R1.size();++i)
682  {
683  R1[i].reinit(block_info.local);
684  R2[i].reinit(block_info.local);
685  }
686  }
687 
688 
689  template <typename number>
690  inline void
691  ResidualLocalBlocksToGlobalBlocks<number>::assemble(
692  FEVector& global,
693  const BlockVector<number>& local,
694  const std::vector<unsigned int>& dof)
695  {
696  for (unsigned int b=0;b<local.n_blocks();++b)
697  for (unsigned int j=0;j<local.block(b).size();++j)
698  {
699  // The coordinates of
700  // the current entry in
701  // DoFHandler
702  // numbering, which
703  // differs from the
704  // block-wise local
705  // numbering we use in
706  // our local vectors
707  const unsigned int jcell = this->block_info.local.local_to_global(b, j);
708  global(dof[jcell]) += local.block(b)(j);
709  }
710  }
711 
712 
713  template <typename number>
714  template <int dim>
715  inline void
716  ResidualLocalBlocksToGlobalBlocks<number>::assemble(
717  const InfoObjects::DoFInfo<dim>& info)
718  {
719  for (unsigned int i=0;i<residuals.n_vectors();++i)
720  {
721  assemble(residuals.vector(i), R1[i], info.indices);
722  R1[i] = 0.;
723  }
724  }
725 
726 
727  template <typename number>
728  template <int dim>
729  inline void
730  ResidualLocalBlocksToGlobalBlocks<number>::assemble(
731  const InfoObjects::DoFInfo<dim>& info1,
732  const InfoObjects::DoFInfo<dim>& info2)
733  {
734  for (unsigned int i=0;i<residuals.n_vectors();++i)
735  {
736  assemble(residuals.vector(i), R1[i], info1.indices);
737  assemble(residuals.vector(i), R2[i], info2.indices);
738  R1[i] = 0.;
739  R2[i] = 0.;
740  }
741  }
742 
743 
744 //----------------------------------------------------------------------//
745 
746 
747  template <typename number>
748  inline void
749  LocalMatrixBlocks<number>::initialize_local(
750  MatrixBlock<FullMatrix<number> >& M,
751  unsigned int row, unsigned int col,
752  unsigned int n_rows, unsigned int n_cols)
753  {
754  M.row = row;
755  M.column = col;
756  M.matrix.reinit(n_rows, n_cols);
757  }
758 
759  template <typename number>
760  template <class MatrixPtr>
761  inline void
762  LocalMatrixBlocks<number>::initialize(
763  const BlockInfo& bi,
764  std::vector<MatrixPtr>& matrices)
765  {
766  block_info = bi;
767 
768  M11.resize(matrices.size());
769  M12.resize(matrices.size());
770  M21.resize(matrices.size());
771  M22.resize(matrices.size());
772 
773  for (unsigned int i=0;i<matrices.size();++i)
774  {
775  const unsigned int row = matrices[i]->row;
776  const unsigned int col = matrices[i]->column;
777  const unsigned int n_rows = block_info.local.block_size(row);
778  const unsigned int n_cols = block_info.local.block_size(col);
779 
780  initialize_local(M11[i], row, col, n_rows, n_cols);
781  initialize_local(M12[i], row, col, n_rows, n_cols);
782  initialize_local(M21[i], row, col, n_rows, n_cols);
783  initialize_local(M22[i], row, col, n_rows, n_cols);
784  }
785  }
786 
787 
788 // ----------------------------------------------------------------------//
789 
790  template <class MATRIX, typename number>
791  inline
792  MatrixLocalBlocksToGlobalBlocks<MATRIX, number>::MatrixLocalBlocksToGlobalBlocks(
793  double threshold)
794  :
795  threshold(threshold)
796  {}
797 
798 
799  template <class MATRIX, typename number>
800  inline void
801  MatrixLocalBlocksToGlobalBlocks<MATRIX, number>::initialize(
802  const BlockInfo& bi,
803  std::vector<MatrixPtr>& m)
804  {
805  LocalMatrixBlocks<number>::initialize(bi, m);
806  matrices = m;
807  }
808 
809 
810 
811  template <class MATRIX, typename number>
812  inline void
813  MatrixLocalBlocksToGlobalBlocks<MATRIX, number>::assemble(
814  MATRIX& global,
815  const FullMatrix<number>& local,
816  unsigned int block_row,
817  unsigned int block_col,
818  const std::vector<unsigned int>& dof1,
819  const std::vector<unsigned int>& dof2)
820  {
821  for (unsigned int j=0;j<local.n_rows();++j)
822  for (unsigned int k=0;k<local.n_cols();++k)
823  if (std::fabs(local(j,k)) >= threshold)
824  {
825  // The coordinates of
826  // the current entry in
827  // DoFHandler
828  // numbering, which
829  // differs from the
830  // block-wise local
831  // numbering we use in
832  // our local matrices
833  const unsigned int jcell = this->block_info.local.local_to_global(block_row, j);
834  const unsigned int kcell = this->block_info.local.local_to_global(block_col, k);
835 
836  // The global dof
837  // indices to assemble
838  // in. Since we may
839  // have face matrices
840  // coupling two
841  // different cells, we
842  // provide two sets of
843  // dof indices.
844  const unsigned int jglobal = this->block_info.global.global_to_local(dof1[jcell]).second;
845  const unsigned int kglobal = this->block_info.global.global_to_local(dof2[kcell]).second;
846 
847  global.add(jglobal, kglobal, local(j,k));
848  }
849  }
850 
851 
852  template <class MATRIX, typename number>
853  template <int dim>
854  inline void
855  MatrixLocalBlocksToGlobalBlocks<MATRIX, number>::assemble(
856  const InfoObjects::DoFInfo<dim>& info)
857  {
858  for (unsigned int i=0;i<matrices.size();++i)
859  {
860  // Row and column index of
861  // the block we are dealing with
862  const unsigned int row = matrices[i]->row;
863  const unsigned int col = matrices[i]->column;
864 
865  assemble(matrices[i]->matrix, this->M11[i].matrix, row, col, info.indices, info.indices);
866  this->M11[i].matrix = 0.;
867  }
868  }
869 
870 
871  template <class MATRIX, typename number>
872  template <int dim>
873  inline void
874  MatrixLocalBlocksToGlobalBlocks<MATRIX, number>::assemble(
875  const InfoObjects::DoFInfo<dim>& info1,
876  const InfoObjects::DoFInfo<dim>& info2)
877  {
878  for (unsigned int i=0;i<matrices.size();++i)
879  {
880  // Row and column index of
881  // the block we are dealing with
882  const unsigned int row = matrices[i]->row;
883  const unsigned int col = matrices[i]->column;
884 
885  assemble(matrices[i]->matrix, this->M11[i].matrix, row, col, info1.indices, info1.indices);
886  assemble(matrices[i]->matrix, this->M12[i].matrix, row, col, info1.indices, info2.indices);
887  assemble(matrices[i]->matrix, this->M21[i].matrix, row, col, info2.indices, info1.indices);
888  assemble(matrices[i]->matrix, this->M22[i].matrix, row, col, info2.indices, info2.indices);
889 
890  this->M11[i].matrix = 0.;
891  this->M12[i].matrix = 0.;
892  this->M21[i].matrix = 0.;
893  this->M22[i].matrix = 0.;
894  }
895  }
896 
897 
898 // ----------------------------------------------------------------------//
899 
900  template <class MATRIX, typename number>
901  inline
902  MGMatrixLocalBlocksToGlobalBlocks<MATRIX, number>::MGMatrixLocalBlocksToGlobalBlocks(
903  double threshold)
904  :
905  threshold(threshold)
906  {}
907 
908 
909  template <class MATRIX, typename number>
910  inline void
911  MGMatrixLocalBlocksToGlobalBlocks<MATRIX, number>::initialize(
912  const BlockInfo& bi,
913  std::vector<MatrixPtr>& m)
914  {
915  LocalMatrixBlocks<number>::initialize(bi, m);
916  matrices = m;
917  }
918 
919 
920  template <class MATRIX, typename number>
921  inline void
922  MGMatrixLocalBlocksToGlobalBlocks<MATRIX, number>::initialize_edge_flux(
923  std::vector<MatrixPtr>& up,
924  std::vector<MatrixPtr>& down)
925  {
926  flux_up = up;
927  flux_down = down;
928  }
929 
930 
931  template <class MATRIX, typename number>
932  inline void
933  MGMatrixLocalBlocksToGlobalBlocks<MATRIX, number>::assemble(
934  MATRIX& global,
935  const FullMatrix<number>& local,
936  unsigned int block_row,
937  unsigned int block_col,
938  const std::vector<unsigned int>& dof1,
939  const std::vector<unsigned int>& dof2,
940  unsigned int level1,
941  unsigned int level2,
942  bool transpose)
943  {
944  for (unsigned int j=0;j<local.n_rows();++j)
945  for (unsigned int k=0;k<local.n_cols();++k)
946  if (std::fabs(local(j,k)) >= threshold)
947  {
948  // The coordinates of
949  // the current entry in
950  // DoFHandler
951  // numbering, which
952  // differs from the
953  // block-wise local
954  // numbering we use in
955  // our local matrices
956  const unsigned int jcell = this->block_info.local.local_to_global(block_row, j);
957  const unsigned int kcell = this->block_info.local.local_to_global(block_col, k);
958 
959  // The global dof
960  // indices to assemble
961  // in. Since we may
962  // have face matrices
963  // coupling two
964  // different cells, we
965  // provide two sets of
966  // dof indices.
967  const unsigned int jglobal = this->block_info.levels[level1].global_to_local(dof1[jcell]).second;
968  const unsigned int kglobal = this->block_info.levels[level2].global_to_local(dof2[kcell]).second;
969 
970  if (transpose)
971  global.add(kglobal, jglobal, local(j,k));
972  else
973  global.add(jglobal, kglobal, local(j,k));
974  }
975  }
976 
977 
978  template <class MATRIX, typename number>
979  template <int dim>
980  inline void
981  MGMatrixLocalBlocksToGlobalBlocks<MATRIX, number>::assemble(
982  const InfoObjects::DoFInfo<dim>& info)
983  {
984  const unsigned int level = info.cell->level();
985 
986  for (unsigned int i=0;i<matrices.size();++i)
987  {
988  // Row and column index of
989  // the block we are dealing with
990  const unsigned int row = matrices[i]->row;
991  const unsigned int col = matrices[i]->column;
992 
993  assemble(matrices[i]->matrix[level], this->M11[i].matrix, row, col,
994  info.indices, info.indices, level, level);
995 
996  this->M11[i].matrix = 0.;
997  }
998  }
999 
1000 
1001  template <class MATRIX, typename number>
1002  template <int dim>
1003  inline void
1004  MGMatrixLocalBlocksToGlobalBlocks<MATRIX, number>::assemble(
1005  const InfoObjects::DoFInfo<dim>& info1,
1006  const InfoObjects::DoFInfo<dim>& info2)
1007  {
1008  const unsigned int level1 = info1.cell->level();
1009  const unsigned int level2 = info2.cell->level();
1010 
1011  for (unsigned int i=0;i<matrices.size();++i)
1012  {
1013  // Row and column index of
1014  // the block we are dealing with
1015  const unsigned int row = matrices[i]->row;
1016  const unsigned int col = matrices[i]->column;
1017 
1018  if (level1 == level2)
1019  {
1020  assemble(matrices[i]->matrix[level1], this->M11[i].matrix, row, col, info1.indices, info1.indices, level1, level1);
1021  assemble(matrices[i]->matrix[level1], this->M12[i].matrix, row, col, info1.indices, info2.indices, level1, level2);
1022  assemble(matrices[i]->matrix[level1], this->M21[i].matrix, row, col, info2.indices, info1.indices, level2, level1);
1023  assemble(matrices[i]->matrix[level1], this->M22[i].matrix, row, col, info2.indices, info2.indices, level2, level2);
1024  }
1025  else
1026  {
1027  Assert(level1 > level2, ExcNotImplemented());
1028  if (flux_up.size() != 0)
1029  {
1030  // Do not add M22,
1031  // which is done by
1032  // the coarser cell
1033  assemble(matrices[i]->matrix[level1], this->M11[i].matrix, row, col,
1034  info1.indices, info1.indices, level1, level1);
1035  assemble(flux_up[i]->matrix[level1], this->M12[i].matrix, row, col,
1036  info1.indices, info2.indices, level1, level2, true);
1037  assemble(flux_down[i]->matrix[level1], this->M21[i].matrix, row, col,
1038  info2.indices, info1.indices, level2, level1);
1039  }
1040  }
1041 
1042  this->M11[i].matrix = 0.;
1043  this->M12[i].matrix = 0.;
1044  this->M21[i].matrix = 0.;
1045  this->M22[i].matrix = 0.;
1046  }
1047  }
1048  }
1049 }
1050 
1051 
1052 #endif
Generated on Fri Aug 30 2013 17:15:54 for OpenFCST: The open-source Fuel Cell Simulation Toolbox by   doxygen 1.8.2