ug4
elem_disc_assemble_util.h
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2010-2015: G-CSC, Goethe University Frankfurt
3  * Author: Andreas Vogel
4  *
5  * This file is part of UG4.
6  *
7  * UG4 is free software: you can redistribute it and/or modify it under the
8  * terms of the GNU Lesser General Public License version 3 (as published by the
9  * Free Software Foundation) with the following additional attribution
10  * requirements (according to LGPL/GPL v3 §7):
11  *
12  * (1) The following notice must be displayed in the Appropriate Legal Notices
13  * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14  *
15  * (2) The following notice must be displayed at a prominent place in the
16  * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17  *
18  * (3) The following bibliography is recommended for citation and must be
19  * preserved in all covered files:
20  * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21  * parallel geometric multigrid solver on hierarchically distributed grids.
22  * Computing and visualization in science 16, 4 (2013), 151-164"
23  * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24  * flexible software system for simulating pde based models on high performance
25  * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26  *
27  * This program is distributed in the hope that it will be useful,
28  * but WITHOUT ANY WARRANTY; without even the implied warranty of
29  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30  * GNU Lesser General Public License for more details.
31  */
32 
33 #ifndef __H__UG__LIB_DISC__SPATIAL_DISC__ELEM_DISC__ELEM_DISC_ASSEMBLE_UTIL__
34 #define __H__UG__LIB_DISC__SPATIAL_DISC__ELEM_DISC__ELEM_DISC_ASSEMBLE_UTIL__
35 
36 // extern includes
37 #include <iostream>
38 #include <vector>
39 
40 // other ug4 modules
41 #include "common/common.h"
42 
43 // intern headers
44 #include "../../reference_element/reference_element.h"
45 #include "./elem_disc_interface.h"
50 
51 #define PROFILE_ELEM_LOOP
52 #ifdef PROFILE_ELEM_LOOP
53  #define EL_PROFILE_FUNC() PROFILE_FUNC()
54  #define EL_PROFILE_BEGIN(name) PROFILE_BEGIN(name)
55  #define EL_PROFILE_END() PROFILE_END()
56 #else
57  #define EL_PROFILE_FUNC()
58  #define EL_PROFILE_BEGIN(name)
59  #define EL_PROFILE_END()
60 #endif
61 
62 
63 namespace ug {
64 
66 
68 
76 template <typename TDomain, typename TAlgebra>
78 {
80  typedef TDomain domain_type;
81 
83  typedef TAlgebra algebra_type;
84 
86  typedef typename algebra_type::vector_type vector_type;
87 
89  typedef typename algebra_type::matrix_type matrix_type;
90 
92 // Assemble Stiffness Matrix
94 
95 public:
112  template <typename TElem, typename TIterator>
113  static void
114  AssembleStiffnessMatrix( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
117  TIterator iterBegin,
118  TIterator iterEnd,
119  int si, bool bNonRegularGrid,
120  matrix_type& A,
121  const vector_type& u,
123  {
124  // check if there are any elements at all, otherwise return immediately
125  if(iterBegin == iterEnd) return;
126 
127  // reference object id
129 
130  // storage for corner coordinates
132 
133  try
134  {
136  vElemDisc, dd->function_pattern(), bNonRegularGrid);
137 
138  // prepare element loop
139  Eval.prepare_elem_loop(id, si);
140 
141  // local indices and local algebra
142  LocalIndices ind; LocalVector locU; LocalMatrix locA;
143 
144  // Loop over all elements
145  for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
146  {
147  // get Element
148  TElem* elem = *iter;
149 
150  // get corner coordinates
151  FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
152 
153  // check if elem is skipped from assembling
154  if(!spAssTuner->element_used(elem)) continue;
155 
156  // get global indices
157  dd->indices(elem, ind, Eval.use_hanging());
158 
159  // adapt local algebra
160  locU.resize(ind); locA.resize(ind);
161 
162  // read local values of u
163  GetLocalVector(locU, u);
164 
165  // prepare element
166  try
167  {
168  Eval.prepare_elem(locU, elem, id, vCornerCoords, ind, true);
169  }
170  UG_CATCH_THROW("AssembleStiffnessMatrix Cannot prepare element.");
171 
172  // Assemble JA
173  locA = 0.0;
174  try
175  {
176  Eval.add_jac_A_elem(locA, locU, elem, vCornerCoords);
177  }
178  UG_CATCH_THROW("AssembleStiffnessMatrix: Cannot compute Jacobian (A).");
179 
180  // send local to global matrix
181  try{
182  spAssTuner->add_local_mat_to_global(A,locA,dd);
183  }
184  UG_CATCH_THROW("AssembleStiffnessMatrix: Cannot add local matrix.");
185  }
186 
187  // finish element loop
188  try
189  {
190  Eval.finish_elem_loop();
191  }
192  UG_CATCH_THROW("AssembleStiffnessMatrix: Cannot finish element loop.");
193 
194  }
195  UG_CATCH_THROW("AssembleStiffnessMatrix': Cannot create Data Evaluator.");
196  }
197 
199 // Assemble Mass Matrix
201 
202 public:
219  template <typename TElem, typename TIterator>
220  static void
221  AssembleMassMatrix( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
224  TIterator iterBegin,
225  TIterator iterEnd,
226  int si, bool bNonRegularGrid,
227  matrix_type& M,
228  const vector_type& u,
230  {
231  // check if there are any elements at all, otherwise return immediately
232  if(iterBegin == iterEnd) return;
233 
234  // reference object id
236 
237  // storage for corner coordinates
239 
240  // prepare for given elem discs
241  try
242  {
244  vElemDisc, dd->function_pattern(), bNonRegularGrid);
245 
246  // prepare element loop
247  Eval.prepare_elem_loop(id, si);
248 
249  // local indices and local algebra
250  LocalIndices ind; LocalVector locU; LocalMatrix locM;
251 
252  // Loop over all elements
253  for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
254  {
255  // get Element
256  TElem* elem = *iter;
257 
258  // get corner coordinates
259  FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
260 
261  // check if elem is skipped from assembling
262  if(!spAssTuner->element_used(elem)) continue;
263 
264  // get global indices
265  dd->indices(elem, ind, Eval.use_hanging());
266 
267  // adapt local algebra
268  locU.resize(ind); locM.resize(ind);
269 
270  // read local values of u
271  GetLocalVector(locU, u);
272 
273  // prepare element
274  try
275  {
276  Eval.prepare_elem(locU, elem, id, vCornerCoords, ind, true);
277  }
278  UG_CATCH_THROW("AssembleMassMatrix: Cannot prepare element.");
279 
280  // Assemble JM
281  locM = 0.0;
282  try
283  {
284  Eval.add_jac_M_elem(locM, locU, elem, vCornerCoords);
285  }
286  UG_CATCH_THROW("AssembleMassMatrix: Cannot compute Jacobian (M).");
287 
288  // send local to global matrix
289  try{
290  spAssTuner->add_local_mat_to_global(M, locM, dd);
291  }
292  UG_CATCH_THROW("AssembleMassMatrix: Cannot add local matrix.");
293  }
294 
295  // finish element loop
296  try
297  {
298  Eval.finish_elem_loop();
299  }
300  UG_CATCH_THROW("AssembleMassMatrix: Cannot finish element loop.");
301 
302  }
303  UG_CATCH_THROW("AssembleMassMatrix: Cannot create Data Evaluator.");
304  }
305 
307 // Assemble (stationary) Jacobian
309 
310 public:
326  template <typename TElem, typename TIterator>
327  static void
328  AssembleJacobian( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
331  TIterator iterBegin,
332  TIterator iterEnd,
333  int si, bool bNonRegularGrid,
334  matrix_type& J,
335  const vector_type& u,
337  {
338  // check if there are any elements at all, otherwise return immediately
339  if(iterBegin == iterEnd) return;
340 
341  // reference object id
343 
344  // storage for corner coordinates
346 
347  // prepare for given elem discs
348  try
349  {
351  vElemDisc, dd->function_pattern(), bNonRegularGrid);
352 
353  // prepare element loop
354  Eval.prepare_elem_loop(id, si);
355 
356  // local indices and local algebra
357  LocalIndices ind; LocalVector locU; LocalMatrix locJ;
358 
359  // Loop over all elements
360  for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
361  {
362  // get Element
363  TElem* elem = *iter;
364 
365  // get corner coordinates
366  FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
367 
368  // check if elem is skipped from assembling
369  if(!spAssTuner->element_used(elem)) continue;
370 
371  // get global indices
372  dd->indices(elem, ind, Eval.use_hanging());
373 
374  // adapt local algebra
375  locU.resize(ind); locJ.resize(ind);
376 
377  // read local values of u
378  GetLocalVector(locU, u);
379 
380  // prepare element
381  try
382  {
383  Eval.prepare_elem(locU, elem, id, vCornerCoords, ind, true);
384  }
385  UG_CATCH_THROW("(stationary) AssembleJacobian: Cannot prepare element.");
386 
387  // reset local algebra
388  locJ = 0.0;
389 
390  // Assemble JA
391  try
392  {
393  Eval.add_jac_A_elem(locJ, locU, elem, vCornerCoords);
394  }
395  UG_CATCH_THROW("(stationary) AssembleJacobian: Cannot compute Jacobian (A).");
396 
397  // send local to global matrix
398  try{
399  spAssTuner->add_local_mat_to_global(J, locJ, dd);
400  }
401  UG_CATCH_THROW("(stationary) AssembleJacobian: Cannot add local matrix.");
402  }
403 
404  // finish element loop
405  try
406  {
407  Eval.finish_elem_loop();
408  }
409  UG_CATCH_THROW("(stationary) AssembleJacobian: Cannot finish element loop.");
410 
411  }
412  UG_CATCH_THROW("(stationary) AssembleJacobian: Cannot create Data Evaluator.");
413  }
414 
416 // Assemble (instationary) Jacobian
418 
419 public:
438  template <typename TElem, typename TIterator>
439  static void
440  AssembleJacobian( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
443  TIterator iterBegin,
444  TIterator iterEnd,
445  int si, bool bNonRegularGrid,
446  matrix_type& J,
448  number s_a0,
450  {
451  // check if there are any elements at all, otherwise return immediately
452  if(iterBegin == iterEnd) return;
453 
454  // reference object id
456 
457  // storage for corner coordinates
459 
460  // get current time and vector
461  const vector_type& u = *vSol->solution(0);
462 
463  // create local time series
464  LocalVectorTimeSeries locTimeSeries;
465  locTimeSeries.read_times(vSol);
466 
467  // prepare for given elem discs
468  try
469  {
471  vElemDisc, dd->function_pattern(), bNonRegularGrid,
472  &locTimeSeries);
473  Eval.set_time_point(0);
474 
475  // prepare element loop
476  Eval.prepare_elem_loop(id, si);
477 
478  // local algebra
479  LocalIndices ind; LocalVector locU; LocalMatrix locJ;
480 
481  EL_PROFILE_BEGIN(Elem_AssembleJacobian);
482  // Loop over all elements
483  for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
484  {
485  // get Element
486  TElem* elem = *iter;
487 
488  // get corner coordinates
489  FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
490 
491  // check if elem is skipped from assembling
492  if(!spAssTuner->element_used(elem)) continue;
493 
494  // get global indices
495  dd->indices(elem, ind, Eval.use_hanging());
496 
497  // adapt local algebra
498  locU.resize(ind); locJ.resize(ind);
499 
500  // read local values of u
501  GetLocalVector(locU, u);
502 
503  // read local values of time series
504  if(Eval.time_series_needed())
505  locTimeSeries.read_values(vSol, ind);
506 
507  // prepare element
508  try
509  {
510  Eval.prepare_elem(locU, elem, id, vCornerCoords, ind, true);
511  }
512  UG_CATCH_THROW("(instationary) AssembleJacobian: Cannot prepare element.");
513 
514  // reset local algebra
515  locJ = 0.0;
516 
517  //EL_PROFILE_BEGIN(Elem_add_JA);
518  // Assemble JA
519  try
520  {
521  Eval.add_jac_A_elem(locJ, locU, elem, vCornerCoords, PT_INSTATIONARY);
522  locJ *= s_a0;
523 
524  Eval.add_jac_A_elem(locJ, locU, elem, vCornerCoords, PT_STATIONARY);
525  }
526  UG_CATCH_THROW("(instationary) AssembleJacobian: Cannot compute Jacobian (A).");
527  //EL_PROFILE_END();
528 
529  // Assemble JM
530  try
531  {
532  Eval.add_jac_M_elem(locJ, locU, elem, vCornerCoords, PT_INSTATIONARY);
533  }
534  UG_CATCH_THROW("(instationary) AssembleJacobian: Cannot compute Jacobian (M).");
535 
536  // send local to global matrix
537  try{
538  spAssTuner->add_local_mat_to_global(J, locJ, dd);
539  }
540  UG_CATCH_THROW("(instationary) AssembleJacobian: Cannot add local matrix.");
541 
542  }
543  EL_PROFILE_END();
544 
545  // finish element loop
546  try
547  {
548  Eval.finish_elem_loop();
549  }
550  UG_CATCH_THROW("(instationary) AssembleJacobian: Cannot finish element loop.");
551 
552  }
553  UG_CATCH_THROW("(instationary) AssembleJacobian: Cannot create Data Evaluator.");
554  }
555 
557 // Assemble (stationary) Defect
559 
560 public:
577  template <typename TElem, typename TIterator>
578  static void
579  AssembleDefect( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
582  TIterator iterBegin,
583  TIterator iterEnd,
584  int si, bool bNonRegularGrid,
585  vector_type& d,
586  const vector_type& u,
588  {
589  // check if at least one element exists, else return
590  if(iterBegin == iterEnd) return;
591 
592  // reference object id
594 
595  // storage for corner coordinates
597 
598  // prepare for given elem discs
599  try
600  {
602  vElemDisc, dd->function_pattern(), bNonRegularGrid);
603 
604  // prepare element loop
605  Eval.prepare_elem_loop(id, si);
606 
607  // local indices and local algebra
608  LocalIndices ind; LocalVector locU, locD, tmpLocD;
609 
610  // Loop over all elements
611  for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
612  {
613  // get Element
614  TElem* elem = *iter;
615 
616  // get corner coordinates
617  FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
618 
619  // check if elem is skipped from assembling
620  if(!spAssTuner->element_used(elem)) continue;
621 
622  // get global indices
623  dd->indices(elem, ind, Eval.use_hanging());
624 
625  // adapt local algebra
626  locU.resize(ind); locD.resize(ind); tmpLocD.resize(ind);
627 
628  // read local values of u
629  GetLocalVector(locU, u);
630 
631  // prepare element
632  try
633  {
634  Eval.prepare_elem(locU, elem, id, vCornerCoords, ind);
635  }
636  UG_CATCH_THROW("(stationary) AssembleDefect: Cannot prepare element.");
637 
638  // ANALOG to 'domain_disc_elem()' - modifies the solution, used
639  // for computing the defect
640  if( spAssTuner->modify_solution_enabled() )
641  {
642  LocalVector& modLocU = locU;
643  try{
644  spAssTuner->modify_LocalSol(modLocU, locU, dd);
645  } UG_CATCH_THROW("Cannot modify local solution.");
646 
647  // recopy modified LocalVector:
648  locU = modLocU;
649  }
650 
651  // reset local algebra
652  locD = 0.0;
653 
654  // Assemble A
655  try
656  {
657  Eval.add_def_A_elem(locD, locU, elem, vCornerCoords);
658  }
659  UG_CATCH_THROW("(stationary) AssembleDefect: Cannot compute Defect (A).");
660 
661  // Assemble rhs
662  try
663  {
664  tmpLocD = 0.0;
665  Eval.add_rhs_elem(tmpLocD, elem, vCornerCoords);
666  locD.scale_append(-1, tmpLocD);
667 
668  }
669  UG_CATCH_THROW("(stationary) AssembleDefect: Cannot compute Rhs.");
670 
671  // send local to global defect
672  try{
673  spAssTuner->add_local_vec_to_global(d, locD, dd);
674  }
675  UG_CATCH_THROW("(stationary) AssembleDefect: Cannot add local vector.");
676  }
677 
678  // finish element loop
679  try
680  {
681  Eval.finish_elem_loop();
682  }
683  UG_CATCH_THROW("(stationary) AssembleDefect: Cannot finish element loop.");
684 
685  }
686  UG_CATCH_THROW("(stationary) AssembleDefect: Cannot create Data Evaluator.");
687  }
688 
690 // Assemble (instationary) Defect
692 
693 public:
712  template <typename TElem, typename TIterator>
713  static void
714  AssembleDefect( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
717  TIterator iterBegin,
718  TIterator iterEnd,
719  int si, bool bNonRegularGrid,
720  vector_type& d,
722  const std::vector<number>& vScaleMass,
723  const std::vector<number>& vScaleStiff,
725  {
726  // check if there are any elements at all, otherwise return immediately
727  if(iterBegin == iterEnd) return;
728 
729  // reference object id
731 
732  // storage for corner coordinates
734 
735  // check time scheme
736  if(vScaleMass.size() != vScaleStiff.size())
737  UG_THROW("(instationary) AssembleDefect: s_a and s_m must have same size.");
738 
739  if(vSol->size() < vScaleStiff.size())
740  UG_THROW("(instationary) AssembleDefect: Time stepping scheme needs at "
741  "least "<<vScaleStiff.size()<<" time steps, but only "<<
742  vSol->size() << " passed.");
743 
744  // create local time series
745  LocalVectorTimeSeries locTimeSeries;
746  locTimeSeries.read_times(vSol);
747 
748  // prepare for given elem discs
749  try
750  {
752  vElemDisc, dd->function_pattern(), bNonRegularGrid,
753  &locTimeSeries, &vScaleMass, &vScaleStiff);
754 
755  // prepare element loop
756  Eval.prepare_elem_loop(id, si);
757 
758  // local indices and local algebra
759  LocalIndices ind; LocalVector locD, tmpLocD;
760 
761  // Loop over all elements
762  for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
763  {
764  // get Element
765  TElem* elem = *iter;
766 
767  // get corner coordinates
768  FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
769 
770  // check if elem is skipped from assembling
771  if(!spAssTuner->element_used(elem)) continue;
772 
773  // get global indices
774  dd->indices(elem, ind, Eval.use_hanging());
775 
776  // adapt local algebra
777  locD.resize(ind); tmpLocD.resize(ind);
778 
779  // read local values of time series
780  locTimeSeries.read_values(vSol, ind);
781 
782  // reset contribution of this element
783  locD = 0.0;
784 
785  // loop all time points and assemble them
786  for(size_t t = 0; t < vScaleStiff.size(); ++t)
787  {
788  number scale_stiff = vScaleStiff[t];
789 
790  // get local solution at timepoint
791  LocalVector& locU = locTimeSeries.solution(t);
792  Eval.set_time_point(t);
793 
794  // prepare element
795  try
796  {
797  Eval.prepare_elem(locU, elem, id, vCornerCoords, ind, false);
798  }
799  UG_CATCH_THROW("(instationary) AssembleDefect: Cannot prepare element.");
800 
801  // Assemble M
802  try
803  {
804  tmpLocD = 0.0;
805  Eval.add_def_M_elem(tmpLocD, locU, elem, vCornerCoords, PT_INSTATIONARY);
806  locD.scale_append(vScaleMass[t], tmpLocD);
807  }
808  UG_CATCH_THROW("(instationary) AssembleDefect: Cannot compute Defect (M).");
809 
810  // Assemble A
811  try
812  {
813  if(scale_stiff != 0.0)
814  {
815  tmpLocD = 0.0;
816  Eval.add_def_A_elem(tmpLocD, locU, elem, vCornerCoords, PT_INSTATIONARY);
817  locD.scale_append(scale_stiff, tmpLocD);
818  }
819 
820  if(t == 0)
821  Eval.add_def_A_elem(locD, locU, elem, vCornerCoords, PT_STATIONARY);
822  }
823  UG_CATCH_THROW("(instationary) AssembleDefect: Cannot compute Defect (A).");
824 
825 
826  // Assemble defect for explicit reaction_rate, reaction and source
827  if( t == 1 ) // only valid at lowest timediscretization order
828  {
829  tmpLocD = 0.0;
830  try
831  {
832  Eval.add_def_A_expl_elem(tmpLocD, locU, elem, vCornerCoords, PT_INSTATIONARY);
833  }
834  UG_CATCH_THROW("(instationary) AssembleDefect explizit: Cannot compute Defect (A).");
835 
836 // UG_ASSERT(vScaleStiff.size() == 2, "Only one step method supported.");
837  const number dt = vSol->time(0)-vSol->time(1);
838  locD.scale_append(dt, tmpLocD);
839  }
840 
841 
842  // Assemble rhs
843  try
844  {
845  if(scale_stiff != 0.0)
846  {
847  tmpLocD = 0.0;
848  Eval.add_rhs_elem(tmpLocD, elem, vCornerCoords, PT_INSTATIONARY);
849  locD.scale_append( - scale_stiff, tmpLocD);
850  }
851 
852  if(t == 0)
853  {
854  tmpLocD = 0.0;
855  Eval.add_rhs_elem(tmpLocD, elem, vCornerCoords, PT_STATIONARY);
856  locD.scale_append( -1.0, tmpLocD);
857  }
858  }
859  UG_CATCH_THROW("(instationary) AssembleDefect: Cannot compute Rhs.");
860  }
861 
862  // send local to global defect
863  try{
864  spAssTuner->add_local_vec_to_global(d, locD, dd);
865  }
866  UG_CATCH_THROW("(instationary) AssembleDefect: Cannot add local vector.");
867  }
868 
869  // finish element loop
870  try
871  {
872  Eval.finish_elem_loop();
873  }
874  UG_CATCH_THROW("(instationary) AssembleDefect: Cannot finish element loop.");
875 
876  }
877  UG_CATCH_THROW("(instationary) AssembleDefect: Cannot create Data Evaluator.");
878  }
879 
881 // Assemble (stationary) Linear problem
883 
884 public:
902  template <typename TElem, typename TIterator>
903  static void
904  AssembleLinear( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
907  TIterator iterBegin,
908  TIterator iterEnd,
909  int si, bool bNonRegularGrid,
910  matrix_type& A,
911  vector_type& rhs,
913  {
914  // check if there are any elements at all, otherwise return immediately
915  if(iterBegin == iterEnd) return;
916 
917  // reference object id
919 
920  // storage for corner coordinates
922 
923  // prepare for given elem discs
924  try
925  {
927  vElemDisc, dd->function_pattern(), bNonRegularGrid);
928 
929  UG_DLOG(DID_ELEM_DISC_ASSEMBLE_UTIL, 2, ">>OCT_DISC_DEBUG: " << "elem_disc_assemble_util.h: " << "AssembleLinear(): DataEvaluator(): " << id << std::endl);
930 
931  // prepare loop
932  Eval.prepare_elem_loop(id, si);
933 
934  UG_DLOG(DID_ELEM_DISC_ASSEMBLE_UTIL, 2, ">>OCT_DISC_DEBUG: " << "elem_disc_assemble_util.h: " << "AssembleLinear(): prepare_elem_loop(): " << id << std::endl);
935 
936  // local indices and local algebra
937  LocalIndices ind; LocalVector locRhs; LocalMatrix locA;
938 
939  // Loop over all elements
940  for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
941  {
942  // get Element
943  TElem* elem = *iter;
944 
945  // get corner coordinates
946  FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
947 
948  // check if elem is skipped from assembling
949  if(!spAssTuner->element_used(elem)) continue;
950 
951  // get global indices
952  dd->indices(elem, ind, Eval.use_hanging());
953 
954  // adapt local algebra
955  locRhs.resize(ind); locA.resize(ind);
956 
957  // prepare element
958  try
959  {
960  UG_DLOG(DID_ELEM_DISC_ASSEMBLE_UTIL, 2, ">>OCT_DISC_DEBUG: " << "elem_disc_assemble_util.h: " << "AssembleLinear(): prepare_elem(): " << id << std::endl);
961  for(int i = 0; i < 8; ++i)
962  UG_DLOG(DID_ELEM_DISC_ASSEMBLE_UTIL, 2, ">>OCT_DISC_DEBUG: " << "elem_disc_assemble_util.h: " << "AssembleLinear(): prepare_elem(): " << "vCornerCoords " << i << ": " << vCornerCoords[i] << std::endl);
963 
964  Eval.prepare_elem(locRhs, elem, id, vCornerCoords, ind, true);
965  }
966  UG_CATCH_THROW("(stationary) AssembleLinear: Cannot prepare element.");
967 
968  // reset local algebra
969  locA = 0.0;
970  locRhs = 0.0;
971 
972  // Assemble JA
973  try
974  {
975  Eval.add_jac_A_elem(locA, locRhs, elem, vCornerCoords);
976  }
977  UG_CATCH_THROW("(stationary) AssembleLinear: Cannot compute Jacobian (A).");
978 
979  // Assemble rhs
980  try
981  {
982  Eval.add_rhs_elem(locRhs, elem, vCornerCoords);
983  }
984  UG_CATCH_THROW("(stationary) AssembleLinear: Cannot compute Rhs.");
985 
986  // send local to global matrix & rhs
987  try{
988  spAssTuner->add_local_mat_to_global(A, locA, dd);
989  spAssTuner->add_local_vec_to_global(rhs, locRhs, dd);
990  }
991  UG_CATCH_THROW("(stationary) AssembleLinear: Cannot add local vector/matrix.");
992  }
993 
994  // finish element loop
995  try
996  {
997  Eval.finish_elem_loop();
998  }
999  UG_CATCH_THROW("(stationary) AssembleLinear: Cannot finish element loop.");
1000 
1001  }
1002  UG_CATCH_THROW("(stationary) AssembleLinear: Cannot create Data Evaluator.");
1003  }
1004 
1006 // Assemble (instationary) Linear problem
1008 
1009 public:
1030  template <typename TElem, typename TIterator>
1031  static void
1032  AssembleLinear( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
1033  ConstSmartPtr<domain_type> spDomain,
1035  TIterator iterBegin,
1036  TIterator iterEnd,
1037  int si, bool bNonRegularGrid,
1038  matrix_type& A,
1039  vector_type& rhs,
1041  const std::vector<number>& vScaleMass,
1042  const std::vector<number>& vScaleStiff,
1044  {
1045  // check if there are any elements at all, otherwise return immediately
1046  if(iterBegin == iterEnd) return;
1047 
1048  // reference object id
1050 
1051  // storage for corner coordinates
1053 
1054  // check time scheme
1055  if(vScaleMass.size() != vScaleStiff.size())
1056  UG_THROW("(instationary) AssembleLinear: s_a and s_m must have same size.");
1057 
1058  if(vSol->size() < vScaleStiff.size())
1059  UG_THROW("(instationary) AssembleLinear: Time stepping scheme needs at "
1060  "least "<<vScaleStiff.size()<<" time steps, but only "<<
1061  vSol->size() << " passed.");
1062 
1063  // create local time solution
1064  LocalVectorTimeSeries locTimeSeries;
1065  locTimeSeries.read_times(vSol);
1066 
1067  // prepare for given elem discs
1068  try
1069  {
1071  vElemDisc, dd->function_pattern(), bNonRegularGrid,
1072  &locTimeSeries, &vScaleMass, &vScaleStiff);
1073 
1074  // prepare loop
1075  Eval.prepare_elem_loop(id, si);
1076 
1077  // local algebra
1078  LocalIndices ind; LocalVector locRhs, tmpLocRhs; LocalMatrix locA, tmpLocA;
1079 
1080  // Loop over all elements
1081  for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
1082  {
1083  // get Element
1084  TElem* elem = *iter;
1085 
1086  // get corner coordinates
1087  FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
1088 
1089  // check if elem is skipped from assembling
1090  if(!spAssTuner->element_used(elem)) continue;
1091 
1092  // get global indices
1093  dd->indices(elem, ind, Eval.use_hanging());
1094 
1095  // adapt local algebra
1096  locRhs.resize(ind); tmpLocRhs.resize(ind);
1097  locA.resize(ind); tmpLocA.resize(ind);
1098 
1099  // read local values of time series
1100  locTimeSeries.read_values(vSol, ind);
1101  Eval.set_time_point(0);
1102 
1103  // reset element contribution
1104  locA = 0.0; locRhs = 0.0;
1105 
1107  // current time step
1108 
1109  // get local solution at time point
1110  LocalVector& locU = locTimeSeries.solution(0);
1111 
1112  // prepare element
1113  try
1114  {
1115  Eval.prepare_elem(locU, elem, id, vCornerCoords, ind, true);
1116  }
1117  UG_CATCH_THROW("(instationary) AssembleLinear: Cannot prepare element.");
1118 
1119  if (!spAssTuner->matrix_is_const())
1120  {
1121  // Assemble JM
1122  try
1123  {
1124  tmpLocA = 0.0;
1125  Eval.add_jac_M_elem(tmpLocA, locU, elem, vCornerCoords, PT_INSTATIONARY);
1126  locA.scale_append(vScaleMass[0], tmpLocA);
1127  }
1128  UG_CATCH_THROW("(instationary) AssembleLinear: Cannot compute Jacobian (M).");
1129 
1130  // Assemble JA
1131  try
1132  {
1133  if (vScaleStiff[0] != 0.0)
1134  {
1135  tmpLocA = 0.0;
1136  Eval.add_jac_A_elem(tmpLocA, locU, elem, vCornerCoords, PT_INSTATIONARY);
1137  locA.scale_append(vScaleStiff[0], tmpLocA);
1138  }
1139  Eval.add_jac_A_elem(locA, locU, elem, vCornerCoords, PT_STATIONARY);
1140  }
1141  UG_CATCH_THROW("(instationary) AssembleLinear: Cannot compute Jacobian (A).");
1142  }
1143  // Assemble rhs
1144  try
1145  {
1146  if (vScaleStiff[0] != 0.0)
1147  {
1148  tmpLocRhs = 0.0;
1149  Eval.add_rhs_elem(tmpLocRhs, elem, vCornerCoords, PT_INSTATIONARY);
1150  locRhs.scale_append(vScaleStiff[0], tmpLocRhs);
1151  }
1152  Eval.add_rhs_elem(locRhs, elem, vCornerCoords, PT_STATIONARY);
1153  }
1154  UG_CATCH_THROW("(instationary) AssembleLinear: Cannot compute Rhs.");
1155 
1157  // old time steps
1158 
1159  // loop all old time points
1160  for(size_t t = 1; t < vScaleStiff.size(); ++t)
1161  {
1162  // get local solution at time point
1163  LocalVector& locU = locTimeSeries.solution(t);
1164  Eval.set_time_point(t);
1165  number scaleStiff = vScaleStiff[t];
1166 
1167  // prepare element
1168  try
1169  {
1170  Eval.prepare_elem(locU, elem, id, vCornerCoords, ind, false);
1171  }
1172  UG_CATCH_THROW("(instationary) AssembleLinear: Cannot prepare element.");
1173 
1174  // Assemble dM
1175  try
1176  {
1177  tmpLocRhs = 0.0;
1178  Eval.add_def_M_elem(tmpLocRhs, locU, elem, vCornerCoords, PT_INSTATIONARY);
1179  locRhs.scale_append(-vScaleMass[t], tmpLocRhs);
1180  }
1181  UG_CATCH_THROW("(instationary) AssembleLinear: Cannot compute Jacobian (M).");
1182 
1183  // Assemble dA
1184  try
1185  {
1186  if (scaleStiff != 0.0)
1187  {
1188  tmpLocRhs = 0.0;
1189  Eval.add_def_A_elem(tmpLocRhs, locU, elem, vCornerCoords, PT_INSTATIONARY);
1190  locRhs.scale_append(-scaleStiff, tmpLocRhs);
1191  }
1192  }
1193  UG_CATCH_THROW("(instationary) AssembleLinear: Cannot compute Jacobian (A).");
1194 
1195  // Assemble rhs
1196  try
1197  {
1198  if (scaleStiff != 0.0)
1199  {
1200  tmpLocRhs = 0.0;
1201  Eval.add_rhs_elem(tmpLocRhs, elem, vCornerCoords, PT_INSTATIONARY);
1202  locRhs.scale_append(scaleStiff, tmpLocRhs);
1203  }
1204  }
1205  UG_CATCH_THROW("(instationary) AssembleLinear: Cannot compute Rhs.");
1206  }
1207 
1208  // send local to global matrix & rhs
1209  try{
1210  if (!spAssTuner->matrix_is_const())
1211  spAssTuner->add_local_mat_to_global(A, locA, dd);
1212  spAssTuner->add_local_vec_to_global(rhs, locRhs, dd);
1213  }
1214  UG_CATCH_THROW("(instationary) AssembleLinear: Cannot add local vector/matrix.");
1215  }
1216 
1217  // finish element loop
1218  try
1219  {
1220  Eval.finish_elem_loop();
1221  }
1222  UG_CATCH_THROW("(instationary) AssembleLinear: Cannot finish element loop.");
1223 
1224  }
1225  UG_CATCH_THROW("(instationary) AssembleLinear: Cannot create Data Evaluator.");
1226  }
1227 
1229 // Assemble Rhs (of a stationary problem)
1231 
1232 public:
1249  template <typename TElem, typename TIterator>
1250  static void
1251  AssembleRhs( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
1252  ConstSmartPtr<domain_type> spDomain,
1254  TIterator iterBegin,
1255  TIterator iterEnd,
1256  int si, bool bNonRegularGrid,
1257  vector_type& rhs,
1258  const vector_type& u,
1260  {
1261  // check if there are any elements at all, otherwise return immediately
1262  if(iterBegin == iterEnd) return;
1263 
1264  // reference object id
1266 
1267  // storage for corner coordinates
1269 
1270  // prepare for given elem discs
1271  try
1272  {
1274  vElemDisc, dd->function_pattern(), bNonRegularGrid);
1275 
1276  // prepare loop
1277  Eval.prepare_elem_loop(id, si);
1278 
1279  // local indices and local algebra
1280  LocalIndices ind; LocalVector locU, locRhs;
1281 
1282  // Loop over all elements
1283  for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
1284  {
1285  // get Element
1286  TElem* elem = *iter;
1287 
1288  // get corner coordinates
1289  FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
1290 
1291  // check if elem is skipped from assembling
1292  if(!spAssTuner->element_used(elem)) continue;
1293 
1294  // get global indices
1295  dd->indices(elem, ind, Eval.use_hanging());
1296 
1297  // adapt local algebra
1298  locU.resize(ind); locRhs.resize(ind);
1299 
1300  // read local values of u
1301  GetLocalVector(locU, u);
1302 
1303  // prepare element
1304  try
1305  {
1306  Eval.prepare_elem(locU, elem, id, vCornerCoords, ind);
1307  }
1308  UG_CATCH_THROW("AssembleRhs: Cannot prepare element.");
1309 
1310  // reset local algebra
1311  locRhs = 0.0;
1312 
1313  // Assemble rhs
1314  try
1315  {
1316  Eval.add_rhs_elem(locRhs, elem, vCornerCoords);
1317  }
1318  UG_CATCH_THROW("AssembleRhs: Cannot compute Rhs.");
1319 
1320  // send local to global rhs
1321  try{
1322  spAssTuner->add_local_vec_to_global(rhs, locRhs, dd);
1323  }
1324  UG_CATCH_THROW("AssembleRhs: Cannot add local vector.");
1325  }
1326 
1327  // finish element loop
1328  try
1329  {
1330  Eval.finish_elem_loop();
1331  }
1332  UG_CATCH_THROW("AssembleRhs: Cannot finish element loop.");
1333 
1334  }
1335  UG_CATCH_THROW("AssembleRhs: Cannot create Data Evaluator.");
1336  }
1337 
1339 // Assemble (instationary) Rhs
1341 
1342 public:
1361  template <typename TElem, typename TIterator>
1362  static void
1363  AssembleRhs( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
1364  ConstSmartPtr<domain_type> spDomain,
1366  TIterator iterBegin,
1367  TIterator iterEnd,
1368  int si, bool bNonRegularGrid,
1369  vector_type& rhs,
1371  const std::vector<number>& vScaleMass,
1372  const std::vector<number>& vScaleStiff,
1374  {
1375  // check if there are any elements at all, otherwise return immediately
1376  if(iterBegin == iterEnd) return;
1377 
1378  // reference object id
1380 
1381  // storage for corner coordinates
1383 
1384  // check time scheme
1385  if(vScaleMass.size() != vScaleStiff.size())
1386  UG_THROW("(instationary) AssembleRhs: s_a and s_m must have same size.");
1387 
1388  if(vSol->size() < vScaleStiff.size())
1389  UG_THROW("(instationary) AssembleRhs: Time stepping scheme needs at "
1390  "least "<<vScaleStiff.size()<<" time steps, but only "<<
1391  vSol->size() << " passed.");
1392 
1393  // get current time
1394  LocalVectorTimeSeries locTimeSeries;
1395  locTimeSeries.read_times(vSol);
1396 
1397  // prepare for given elem discs
1398  try
1399  {
1401  vElemDisc, dd->function_pattern(), bNonRegularGrid,
1402  &locTimeSeries, &vScaleMass, &vScaleStiff);
1403 
1404  // prepare loop
1405  Eval.prepare_elem_loop(id, si);
1406 
1407  // local algebra
1408  LocalIndices ind; LocalVector locRhs, tmpLocRhs;
1409 
1410  // Loop over all elements
1411  for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
1412  {
1413  // get Element
1414  TElem* elem = *iter;
1415 
1416  // get corner coordinates
1417  FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
1418 
1419  // check if elem is skipped from assembling
1420  if(!spAssTuner->element_used(elem)) continue;
1421 
1422  // get global indices
1423  dd->indices(elem, ind, Eval.use_hanging());
1424 
1425  // adapt local algebra
1426  locRhs.resize(ind); tmpLocRhs.resize(ind);
1427 
1428  // read local values of time series
1429  locTimeSeries.read_values(vSol, ind);
1430  Eval.set_time_point(0);
1431 
1432  // reset element contribution
1433  locRhs = 0.0;
1434 
1436  // current time step
1437 
1438  // get local solution at time point
1439  LocalVector& locU = locTimeSeries.solution(0);
1440 
1441  // prepare element
1442  try
1443  {
1444  Eval.prepare_elem(locU, elem, id, vCornerCoords, ind, false);
1445  }
1446  UG_CATCH_THROW("(instationary) AssembleRhs: Cannot prepare element.");
1447 
1448  // Assemble rhs
1449  try
1450  {
1451  tmpLocRhs = 0.0;
1452  Eval.add_rhs_elem(tmpLocRhs, elem, vCornerCoords, PT_INSTATIONARY);
1453  locRhs.scale_append(vScaleStiff[0], tmpLocRhs);
1454 
1455  Eval.add_rhs_elem(locRhs, elem, vCornerCoords, PT_STATIONARY);
1456  }
1457  UG_CATCH_THROW("(instationary) AssembleRhs: Cannot compute Rhs.");
1458 
1460  // old time steps
1461 
1462  // loop all old time points
1463  for(size_t t = 1; t < vScaleStiff.size(); ++t)
1464  {
1465  // get local solution at time point
1466  LocalVector& locU = locTimeSeries.solution(t);
1467  Eval.set_time_point(t);
1468 
1469  // prepare element
1470  try
1471  {
1472  Eval.prepare_elem(locU, elem, id, vCornerCoords, ind, false);
1473  }
1474  UG_CATCH_THROW("(instationary) AssembleRhs: Cannot prepare element.");
1475 
1476  // Assemble dM
1477  try
1478  {
1479  tmpLocRhs = 0.0;
1480  Eval.add_def_M_elem(tmpLocRhs, locU, elem, vCornerCoords, PT_INSTATIONARY);
1481  locRhs.scale_append(-vScaleMass[t], tmpLocRhs);
1482  }
1483  UG_CATCH_THROW("(instationary) AssembleRhs: Cannot compute Jacobian (M).");
1484 
1485  // Assemble dA
1486  try
1487  {
1488  tmpLocRhs = 0.0;
1489  Eval.add_def_A_elem(tmpLocRhs, locU, elem, vCornerCoords, PT_INSTATIONARY);
1490  locRhs.scale_append(-vScaleStiff[t], tmpLocRhs);
1491  }
1492  UG_CATCH_THROW("(instationary) AssembleRhs: Cannot compute Jacobian (A).");
1493 
1494  // Assemble rhs
1495  try
1496  {
1497  tmpLocRhs = 0.0;
1498  Eval.add_rhs_elem(tmpLocRhs, elem, vCornerCoords, PT_INSTATIONARY);
1499  locRhs.scale_append(vScaleStiff[t], tmpLocRhs);
1500  }
1501  UG_CATCH_THROW("(instationary) AssembleRhs: Cannot compute Rhs.");
1502  }
1503 
1504  // send local to global rhs
1505  try{
1506  spAssTuner->add_local_vec_to_global(rhs, locRhs, dd);
1507  }
1508  UG_CATCH_THROW("(instationary) AssembleRhs: Cannot add local vector.");
1509  }
1510 
1511  // finish element loop
1512  try
1513  {
1514  Eval.finish_elem_loop();
1515  }
1516  UG_CATCH_THROW("(instationary) AssembleRhs: Cannot finish element loop.");
1517 
1518  }
1519  UG_CATCH_THROW("(instationary) AssembleRhs: Cannot create Data Evaluator.");
1520  }
1521 
1522 // //////////////////////////////////////////////////////////////////////////////
1523 // Prepare Timestep (for instationary problems)
1524 // /////////////////////////////////////////////////////////////////////////////
1525 public:
1537  static void
1539  (
1540  const std::vector<IElemDisc<domain_type>*>& vElemDisc,
1542  bool bNonRegularGrid,
1544  number future_time,
1546  )
1547  {
1548  // get current time and vector
1549  const vector_type& u = *vSol->solution(0);
1550 
1551  // create local time series
1552  LocalVectorTimeSeries locTimeSeries;
1553  locTimeSeries.read_times(vSol);
1554 
1555  try
1556  {
1558  vElemDisc, dd->function_pattern(), bNonRegularGrid,
1559  &locTimeSeries);
1560  Eval.set_time_point(0);
1561 
1562  // prepare timestep
1563  try
1564  {
1566  size_t algebra_index = bridge::AlgebraTypeIDProvider::instance().id<algebra_type>();
1567  Eval.prepare_timestep(future_time, vSol->time(0), &vp, algebra_index);
1568  }
1569  UG_CATCH_THROW("(instationary) PrepareTimestep: Cannot prepare timestep.");
1570 
1571  }
1572  UG_CATCH_THROW("(instationary) PrepareTimestep: Cannot create Data Evaluator.");
1573  }
1574 
1576 // Prepare Timestep Elem (for instationary problems)
1578 public:
1595  template <typename TElem, typename TIterator>
1596  static void
1597  PrepareTimestepElem(const std::vector<IElemDisc<domain_type>*>& vElemDisc,
1598  ConstSmartPtr<domain_type> spDomain,
1600  TIterator iterBegin,
1601  TIterator iterEnd,
1602  int si, bool bNonRegularGrid,
1605  {
1606  // check if there are any elements at all, otherwise return immediately
1607  if(iterBegin == iterEnd) return;
1608 
1609  // reference object id
1611 
1612  // storage for corner coordinates
1614 
1615  // get current time and vector
1616  const vector_type& u = *vSol->solution(0);
1617 
1618  // create local time series
1619  LocalVectorTimeSeries locTimeSeries;
1620  locTimeSeries.read_times(vSol);
1621 
1622  try
1623  {
1625  vElemDisc, dd->function_pattern(), bNonRegularGrid,
1626  &locTimeSeries);
1627  Eval.set_time_point(0);
1628 
1629  // prepare element loop
1630  Eval.prepare_elem_loop(id, si);
1631 
1632  // local algebra
1633  LocalIndices ind; LocalVector locU;
1634 
1635  // Loop over all elements
1636  for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
1637  {
1638  // get Element
1639  TElem* elem = *iter;
1640 
1641  // get corner coordinates
1642  FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
1643 
1644  // check if elem is skipped from assembling
1645  if(!spAssTuner->element_used(elem)) continue;
1646 
1647  // get global indices
1648  dd->indices(elem, ind, Eval.use_hanging());
1649 
1650  // adapt local algebra
1651  locU.resize(ind);
1652 
1653  // read local values of u
1654  GetLocalVector(locU, u);
1655 
1656  // read local values of time series
1657  if(Eval.time_series_needed())
1658  locTimeSeries.read_values(vSol, ind);
1659 
1660  // prepare timestep
1661  try
1662  {
1663  Eval.prepare_timestep_elem(vSol->time(0), locU, elem, vCornerCoords);
1664  }
1665  UG_CATCH_THROW("(instationary) PrepareTimestep: Cannot prepare timestep.");
1666  }
1667 
1668  // finish element loop
1669  try
1670  {
1671  Eval.finish_elem_loop();
1672  }
1673  UG_CATCH_THROW("AssembleMassMatrix: Cannot finish element loop.");
1674 
1675  }
1676  UG_CATCH_THROW("(instationary) PrepareTimestep: Cannot create Data Evaluator.");
1677  }
1678 
1679 // //////////////////////////////////////////////////////////////////////////////
1680 // Finish Timestep (for instationary problems)
1681 // /////////////////////////////////////////////////////////////////////////////
1682 public:
1694  static void
1696  (
1697  const std::vector<IElemDisc<domain_type>*>& vElemDisc,
1699  bool bNonRegularGrid,
1702  )
1703  {
1704  // get current time and vector
1705  const vector_type& u = *vSol->solution(0);
1706 
1707  // create local time series
1708  LocalVectorTimeSeries locTimeSeries;
1709  locTimeSeries.read_times(vSol);
1710 
1711  try
1712  {
1714  vElemDisc, dd->function_pattern(), bNonRegularGrid,
1715  &locTimeSeries);
1716  Eval.set_time_point(0);
1717 
1718  // finish time step
1719  try
1720  {
1722  size_t algebra_index = bridge::AlgebraTypeIDProvider::instance().id<algebra_type>();
1723  Eval.finish_timestep(vSol->time(0), &vp, algebra_index);
1724  }
1725  UG_CATCH_THROW("(instationary) FinishTimestep: Cannot prepare time step.");
1726 
1727  }
1728  UG_CATCH_THROW("(instationary) FinishTimestep: Cannot create Data Evaluator.");
1729  }
1730 
1732 // Finish Timestep Elem (for instationary problems)
1734 
1735 public:
1751  template <typename TElem, typename TIterator>
1752  static void
1753  FinishTimestepElem(const std::vector<IElemDisc<domain_type>*>& vElemDisc,
1754  ConstSmartPtr<domain_type> spDomain,
1756  TIterator iterBegin,
1757  TIterator iterEnd,
1758  int si, bool bNonRegularGrid,
1761  {
1762  // check if there are any elements at all, otherwise return immediately
1763  if(iterBegin == iterEnd) return;
1764 
1765  // reference object id
1767 
1768  // storage for corner coordinates
1770 
1771  // get current time and vector
1772  const vector_type& u = *vSol->solution(0);
1773 
1774  // create local time series
1775  LocalVectorTimeSeries locTimeSeries;
1776  locTimeSeries.read_times(vSol);
1777 
1778  // prepare for given elem discs
1779  try
1780  {
1782  vElemDisc, dd->function_pattern(), bNonRegularGrid,
1783  &locTimeSeries);
1784  Eval.set_time_point(0);
1785 
1786  // prepare loop
1787  Eval.prepare_elem_loop(id, si);
1788 
1789  // local algebra
1790  LocalIndices ind; LocalVector locU;
1791 
1792  // Loop over all elements
1793  for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
1794  {
1795  // get Element
1796  TElem* elem = *iter;
1797 
1798  // get corner coordinates
1799  FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
1800 
1801  // check if elem is skipped from assembling
1802  if(!spAssTuner->element_used(elem)) continue;
1803 
1804  // get global indices
1805  dd->indices(elem, ind, Eval.use_hanging());
1806 
1807  // adapt local algebra
1808  locU.resize(ind);
1809 
1810  // read local values of u
1811  GetLocalVector(locU, u);
1812 
1813  // read local values of time series
1814  if(Eval.time_series_needed())
1815  locTimeSeries.read_values(vSol, ind);
1816 
1817  // finish timestep
1818  try
1819  {
1820  Eval.finish_timestep_elem(locTimeSeries.time(0), locU, elem, vCornerCoords);
1821  }
1822  UG_CATCH_THROW("(instationary) FinishTimestepElem: Cannot finish timestep.");
1823  }
1824 
1825  // finish element loop
1826  try
1827  {
1828  Eval.finish_elem_loop();
1829  }
1830  UG_CATCH_THROW("AssembleMassMatrix: Cannot finish element loop.");
1831 
1832  }
1833  UG_CATCH_THROW("(instationary) FinishTimestepElem: Cannot create Data Evaluator");
1834  }
1835 
1837 // Init. all exports (an optional operation, to use the exports for plotting etc.)
1839 
1840 public:
1855  template <typename TElem, typename TIterator>
1856  static void
1857  InitAllExports(const std::vector<IElemDisc<domain_type>*>& vElemDisc,
1859  TIterator iterBegin,
1860  TIterator iterEnd,
1861  int si, bool bNonRegularGrid, bool bAsTimeDependent)
1862  {
1863  // check if there are any elements at all, otherwise return immediately
1864  if(iterBegin == iterEnd) return;
1865 
1866  // reference object id
1868 
1869  // dummy local time series (only to simulate the instationary case for the initialization)
1870  LocalVectorTimeSeries locTimeSeries;
1871 
1872  // prepare for given elem discs
1873  try
1874  {
1876  vElemDisc, dd->function_pattern(), bNonRegularGrid,
1877  bAsTimeDependent? &locTimeSeries : NULL);
1878  Eval.set_time_point(0);
1879 
1880  // prepare loop
1881  Eval.prepare_elem_loop(id, si);
1882 
1883 
1884  // finish element loop
1885  Eval.finish_elem_loop();
1886 
1887  }
1888  UG_CATCH_THROW("InitAllExports: Data Evaluator failed");
1889  }
1890 
1892 // Error estimator (stationary)
1894 
1895 public:
1910  template <typename TElem, typename TIterator>
1911  static void
1913  (
1914  const std::vector<IElemError<domain_type>*>& vElemDisc,
1915  ConstSmartPtr<domain_type> spDomain,
1917  TIterator iterBegin,
1918  TIterator iterEnd,
1919  int si,
1920  bool bNonRegularGrid,
1921  const vector_type& u
1922  )
1923  {
1924  // check if there are any elements at all, otherwise return immediately
1925  if(iterBegin == iterEnd) return;
1926 
1927  // reference object id
1929 
1930  // storage for corner coordinates
1932 
1933  // prepare for given elem discs
1934  try
1935  {
1937  vElemDisc, dd->function_pattern(), bNonRegularGrid);
1938 
1939  // prepare element loop
1940  Eval.prepare_err_est_elem_loop(id, si);
1941 
1942  // local indices and local algebra
1943  LocalIndices ind; LocalVector locU;
1944 
1945  // Loop over all elements
1946  for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
1947  {
1948  // get Element
1949  TElem* elem = *iter;
1950 
1951  // get corner coordinates
1952  FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
1953 
1954  // get global indices
1955  dd->indices(elem, ind, Eval.use_hanging());
1956 
1957  // adapt local algebra
1958  locU.resize(ind);
1959 
1960  // read local values of u
1961  GetLocalVector(locU, u);
1962 
1963  // prepare element
1964  try
1965  {
1966  Eval.prepare_err_est_elem(locU, elem, vCornerCoords, ind, false);
1967  }
1968  UG_CATCH_THROW("(stationary) AssembleRhs: Cannot prepare element.");
1969 
1970  // assemble stiffness part
1971  try
1972  {
1973  Eval.compute_err_est_A_elem(locU, elem, vCornerCoords, ind);
1974  }
1975  UG_CATCH_THROW("AssembleErrorEstimator: Cannot assemble the error estimator for stiffness part.");
1976 
1977  // assemble rhs part
1978  try
1979  {
1980  Eval.compute_err_est_rhs_elem(elem, vCornerCoords, ind);
1981  }
1982  UG_CATCH_THROW("AssembleErrorEstimator: Cannot assemble the error estimator for stiffness part.");
1983 
1984  }
1985 
1986  // finish element loop
1987  try
1988  {
1989  Eval.finish_err_est_elem_loop();
1990  }
1991  UG_CATCH_THROW("AssembleErrorEstimator: Cannot finish element loop.");
1992 
1993  }
1994  UG_CATCH_THROW("AssembleErrorEstimator: Cannot create Data Evaluator.");
1995  }
1996 
1998 // Error estimator (for time-dependent problems)
2000 
2001 public:
2018  template <typename TElem, typename TIterator>
2019  static void
2021  (
2022  const std::vector<IElemError<domain_type>*>& vElemDisc,
2023  ConstSmartPtr<domain_type> spDomain,
2025  TIterator iterBegin,
2026  TIterator iterEnd,
2027  int si,
2028  bool bNonRegularGrid,
2029  const std::vector<number>& vScaleMass,
2030  const std::vector<number>& vScaleStiff,
2032  )
2033  {
2034  // check if there are any elements at all, otherwise return immediately
2035  if (iterBegin == iterEnd) return;
2036 
2037  // reference object id
2039 
2040  // storage for corner coordinates
2042 
2043  // check time scheme
2044  if(vScaleMass.size() != vScaleStiff.size())
2045  UG_THROW("(instationary) AssembleErrorEstimator: s_a and s_m must have same size.");
2046 
2047  if(vSol->size() < vScaleStiff.size())
2048  UG_THROW("(instationary) AssembleErrorEstimator: Time stepping scheme needs at "
2049  "least "<<vScaleStiff.size()<<" time steps, but only "<<
2050  vSol->size() << " passed.");
2051 
2052  // create local time series
2053  LocalVectorTimeSeries locTimeSeries;
2054  locTimeSeries.read_times(vSol);
2055 
2056  // prepare for given elem discs
2057  try
2058  {
2060  vElemDisc, dd->function_pattern(), bNonRegularGrid,
2061  &locTimeSeries, &vScaleMass, &vScaleStiff);
2062 
2063  // prepare element loop
2064  Eval.prepare_err_est_elem_loop(id, si);
2065 
2066  // local indices and local algebra
2067  LocalIndices ind;
2068 
2069  // loop over all elements
2070  for (TIterator iter = iterBegin; iter != iterEnd; ++iter)
2071  {
2072  // get Element
2073  TElem* elem = *iter;
2074 
2075  // get corner coordinates
2076  FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
2077 
2078  // get global indices
2079  dd->indices(elem, ind, Eval.use_hanging());
2080 
2081  // read local values of time series
2082  locTimeSeries.read_values(vSol, ind);
2083 
2084  // loop all time points and assemble them
2085  for (std::size_t t = 0; t < vScaleStiff.size(); ++t)
2086  {
2087  // get local solution at timepoint
2088  LocalVector& locU = locTimeSeries.solution(t);
2089  Eval.set_time_point(t);
2090 
2091  // prepare element
2092  try
2093  {
2094  Eval.prepare_err_est_elem(locU, elem, vCornerCoords, ind, false);
2095  }
2096  UG_CATCH_THROW("AssembleErrorEstimator: Cannot prepare element.");
2097 
2098  // assemble stiffness part
2099  try
2100  {
2101  Eval.compute_err_est_A_elem(locU, elem, vCornerCoords, ind, vScaleMass[t], vScaleStiff[t]);
2102  }
2103  UG_CATCH_THROW("AssembleErrorEstimator: Cannot assemble the error estimator for stiffness part.");
2104 
2105  // assemble mass part
2106  try
2107  {
2108  Eval.compute_err_est_M_elem(locU, elem, vCornerCoords, ind, vScaleMass[t], vScaleStiff[t]);
2109  }
2110  UG_CATCH_THROW("AssembleErrorEstimator: Cannot assemble the error estimator for stiffness part.");
2111 
2112  // assemble rhs part
2113  try
2114  {
2115  Eval.compute_err_est_rhs_elem(elem, vCornerCoords, ind, vScaleMass[t], vScaleStiff[t]);
2116  }
2117  UG_CATCH_THROW("AssembleErrorEstimator: Cannot assemble the error estimator for stiffness part.");
2118  }
2119  }
2120 
2121  // finish element loop
2122  try
2123  {
2124  Eval.finish_err_est_elem_loop();
2125  }
2126  UG_CATCH_THROW("AssembleErrorEstimator: Cannot finish element loop.");
2127 
2128  }
2129  UG_CATCH_THROW("AssembleErrorEstimator: Cannot create Data Evaluator.");
2130  }
2131 
2132 }; // class StdGlobAssembler
2133 
2134 } // end namespace ug
2135 
2136 
2137 #endif /* __H__UG__LIB_DISC__SPATIAL_DISC__ELEM_DISC__ELEM_DISC_ASSEMBLE_UTIL__ */
The AssemblingTuner class combines tools to adapt the assembling routine.
Definition: ass_tuner.h:90
Evaluation for IElemDisc.
Definition: data_evaluator.h:169
Definition: debug_id.h:94
Evaluation for IElemError.
Definition: data_evaluator.h:151
Definition: elem_disc_interface.h:800
Definition: elem_disc_interface.h:760
Definition: local_algebra.h:50
Definition: local_algebra.h:422
void resize(const LocalIndices &ind)
resize for same ind in row and column
Definition: local_algebra.h:445
this_type & scale_append(number s, const this_type &rhs)
add scaled matrix
Definition: local_algebra.h:523
Definition: local_algebra.h:198
this_type & scale_append(number s, const this_type &rhs)
add a scaled vector
Definition: local_algebra.h:292
void resize(const LocalIndices &ind)
resize for current local indices
Definition: local_algebra.h:214
time series of local vectors
Definition: solution_time_series.h:167
void read_values(ConstSmartPtr< VectorTimeSeries< TVector > > vecTimeSeries, LocalIndices &ind)
extract local values from global vectors
Definition: solution_time_series.h:196
void read_times(ConstSmartPtr< VectorTimeSeries< TVector > > vecTimeSeries)
extract time points
Definition: solution_time_series.h:208
number time(size_t i) const
returns time point i
Definition: solution_time_series.h:176
const LocalVector & solution(size_t i) const
returns the local vector for the i'th time point
Definition: solution_time_series.h:182
a mathematical Vector with N entries.
Definition: math_vector.h:97
Global assembler based on the straightforward application of the local discretizations.
Definition: elem_disc_assemble_util.h:78
static void FinishTimestep(const std::vector< IElemDisc< domain_type > * > &vElemDisc, ConstSmartPtr< DoFDistribution > dd, bool bNonRegularGrid, ConstSmartPtr< VectorTimeSeries< vector_type > > vSol, ConstSmartPtr< AssemblingTuner< TAlgebra > > spAssTuner)
Definition: elem_disc_assemble_util.h:1696
static void AssembleMassMatrix(const std::vector< IElemDisc< domain_type > * > &vElemDisc, ConstSmartPtr< domain_type > spDomain, ConstSmartPtr< DoFDistribution > dd, TIterator iterBegin, TIterator iterEnd, int si, bool bNonRegularGrid, matrix_type &M, const vector_type &u, ConstSmartPtr< AssemblingTuner< TAlgebra > > spAssTuner)
Definition: elem_disc_assemble_util.h:221
static void AssembleLinear(const std::vector< IElemDisc< domain_type > * > &vElemDisc, ConstSmartPtr< domain_type > spDomain, ConstSmartPtr< DoFDistribution > dd, TIterator iterBegin, TIterator iterEnd, int si, bool bNonRegularGrid, matrix_type &A, vector_type &rhs, ConstSmartPtr< AssemblingTuner< TAlgebra > > spAssTuner)
Definition: elem_disc_assemble_util.h:904
static void AssembleErrorEstimator(const std::vector< IElemError< domain_type > * > &vElemDisc, ConstSmartPtr< domain_type > spDomain, ConstSmartPtr< DoFDistribution > dd, TIterator iterBegin, TIterator iterEnd, int si, bool bNonRegularGrid, const vector_type &u)
Definition: elem_disc_assemble_util.h:1913
static void AssembleErrorEstimator(const std::vector< IElemError< domain_type > * > &vElemDisc, ConstSmartPtr< domain_type > spDomain, ConstSmartPtr< DoFDistribution > dd, TIterator iterBegin, TIterator iterEnd, int si, bool bNonRegularGrid, const std::vector< number > &vScaleMass, const std::vector< number > &vScaleStiff, ConstSmartPtr< VectorTimeSeries< vector_type > > vSol)
Definition: elem_disc_assemble_util.h:2021
static void AssembleJacobian(const std::vector< IElemDisc< domain_type > * > &vElemDisc, ConstSmartPtr< domain_type > spDomain, ConstSmartPtr< DoFDistribution > dd, TIterator iterBegin, TIterator iterEnd, int si, bool bNonRegularGrid, matrix_type &J, const vector_type &u, ConstSmartPtr< AssemblingTuner< TAlgebra > > spAssTuner)
Definition: elem_disc_assemble_util.h:328
TAlgebra algebra_type
Algebra type.
Definition: elem_disc_assemble_util.h:83
algebra_type::vector_type vector_type
Vector type in the algebra.
Definition: elem_disc_assemble_util.h:86
static void PrepareTimestep(const std::vector< IElemDisc< domain_type > * > &vElemDisc, ConstSmartPtr< DoFDistribution > dd, bool bNonRegularGrid, ConstSmartPtr< VectorTimeSeries< vector_type > > vSol, number future_time, ConstSmartPtr< AssemblingTuner< TAlgebra > > spAssTuner)
Definition: elem_disc_assemble_util.h:1539
static void AssembleRhs(const std::vector< IElemDisc< domain_type > * > &vElemDisc, ConstSmartPtr< domain_type > spDomain, ConstSmartPtr< DoFDistribution > dd, TIterator iterBegin, TIterator iterEnd, int si, bool bNonRegularGrid, vector_type &rhs, const vector_type &u, ConstSmartPtr< AssemblingTuner< TAlgebra > > spAssTuner)
Definition: elem_disc_assemble_util.h:1251
static void PrepareTimestepElem(const std::vector< IElemDisc< domain_type > * > &vElemDisc, ConstSmartPtr< domain_type > spDomain, ConstSmartPtr< DoFDistribution > dd, TIterator iterBegin, TIterator iterEnd, int si, bool bNonRegularGrid, ConstSmartPtr< VectorTimeSeries< vector_type > > vSol, ConstSmartPtr< AssemblingTuner< TAlgebra > > spAssTuner)
Definition: elem_disc_assemble_util.h:1597
static void AssembleDefect(const std::vector< IElemDisc< domain_type > * > &vElemDisc, ConstSmartPtr< domain_type > spDomain, ConstSmartPtr< DoFDistribution > dd, TIterator iterBegin, TIterator iterEnd, int si, bool bNonRegularGrid, vector_type &d, const vector_type &u, ConstSmartPtr< AssemblingTuner< TAlgebra > > spAssTuner)
Definition: elem_disc_assemble_util.h:579
static void AssembleStiffnessMatrix(const std::vector< IElemDisc< domain_type > * > &vElemDisc, ConstSmartPtr< domain_type > spDomain, ConstSmartPtr< DoFDistribution > dd, TIterator iterBegin, TIterator iterEnd, int si, bool bNonRegularGrid, matrix_type &A, const vector_type &u, ConstSmartPtr< AssemblingTuner< TAlgebra > > spAssTuner)
Definition: elem_disc_assemble_util.h:114
static void AssembleLinear(const std::vector< IElemDisc< domain_type > * > &vElemDisc, ConstSmartPtr< domain_type > spDomain, ConstSmartPtr< DoFDistribution > dd, TIterator iterBegin, TIterator iterEnd, int si, bool bNonRegularGrid, matrix_type &A, vector_type &rhs, ConstSmartPtr< VectorTimeSeries< vector_type > > vSol, const std::vector< number > &vScaleMass, const std::vector< number > &vScaleStiff, ConstSmartPtr< AssemblingTuner< TAlgebra > > spAssTuner)
Definition: elem_disc_assemble_util.h:1032
algebra_type::matrix_type matrix_type
Matrix type in the algebra.
Definition: elem_disc_assemble_util.h:89
static void AssembleRhs(const std::vector< IElemDisc< domain_type > * > &vElemDisc, ConstSmartPtr< domain_type > spDomain, ConstSmartPtr< DoFDistribution > dd, TIterator iterBegin, TIterator iterEnd, int si, bool bNonRegularGrid, vector_type &rhs, ConstSmartPtr< VectorTimeSeries< vector_type > > vSol, const std::vector< number > &vScaleMass, const std::vector< number > &vScaleStiff, ConstSmartPtr< AssemblingTuner< TAlgebra > > spAssTuner)
Definition: elem_disc_assemble_util.h:1363
TDomain domain_type
Domain type.
Definition: elem_disc_assemble_util.h:80
static void AssembleJacobian(const std::vector< IElemDisc< domain_type > * > &vElemDisc, ConstSmartPtr< domain_type > spDomain, ConstSmartPtr< DoFDistribution > dd, TIterator iterBegin, TIterator iterEnd, int si, bool bNonRegularGrid, matrix_type &J, ConstSmartPtr< VectorTimeSeries< vector_type > > vSol, number s_a0, ConstSmartPtr< AssemblingTuner< TAlgebra > > spAssTuner)
Definition: elem_disc_assemble_util.h:440
static void AssembleDefect(const std::vector< IElemDisc< domain_type > * > &vElemDisc, ConstSmartPtr< domain_type > spDomain, ConstSmartPtr< DoFDistribution > dd, TIterator iterBegin, TIterator iterEnd, int si, bool bNonRegularGrid, vector_type &d, ConstSmartPtr< VectorTimeSeries< vector_type > > vSol, const std::vector< number > &vScaleMass, const std::vector< number > &vScaleStiff, ConstSmartPtr< AssemblingTuner< TAlgebra > > spAssTuner)
Definition: elem_disc_assemble_util.h:714
static void FinishTimestepElem(const std::vector< IElemDisc< domain_type > * > &vElemDisc, ConstSmartPtr< domain_type > spDomain, ConstSmartPtr< DoFDistribution > dd, TIterator iterBegin, TIterator iterEnd, int si, bool bNonRegularGrid, ConstSmartPtr< VectorTimeSeries< vector_type > > vSol, ConstSmartPtr< AssemblingTuner< TAlgebra > > spAssTuner)
Definition: elem_disc_assemble_util.h:1753
static void InitAllExports(const std::vector< IElemDisc< domain_type > * > &vElemDisc, ConstSmartPtr< DoFDistribution > dd, TIterator iterBegin, TIterator iterEnd, int si, bool bNonRegularGrid, bool bAsTimeDependent)
Definition: elem_disc_assemble_util.h:1857
time series of solutions and corresponding time point
Definition: solution_time_series.h:59
Definition: grid_base_object_traits.h:68
#define EL_PROFILE_BEGIN(name)
Definition: elem_disc_assemble_util.h:54
#define EL_PROFILE_END()
Definition: elem_disc_assemble_util.h:55
#define UG_CATCH_THROW(msg)
Definition: error.h:64
#define UG_THROW(msg)
Definition: error.h:57
#define UG_DLOG(__debugID__, level, msg)
Definition: log.h:298
double number
Definition: types.h:124
SmartPtr< typename TAlgebra::vector_type > Eval(SmartPtr< VecScaleAddClass< TAlgebra > > vsac)
Definition: bridge_mat_vec_operations.h:141
const int NUM_VERTICES
Definition: hexahedron_rules.h:45
the ug namespace
ReferenceObjectID
these ids are used to identify the shape of a geometric object.
Definition: grid_base_objects.h:74
@ PT_INSTATIONARY
Definition: data_evaluator.h:42
@ PT_STATIONARY
Definition: data_evaluator.h:42
@ NONE
Definition: data_import.h:44
@ STIFF
Definition: data_import.h:46
@ RHS
Definition: data_import.h:47
@ EXPL
Definition: data_import.h:48
@ MASS
Definition: data_import.h:45
void FillCornerCoordinates(typename TDomain::position_type vCornerCoordsOut[], const TElem &elem, const TDomain &domain)
returns the corner coordinates of a geometric object
Definition: domain_util_impl.h:114
void GetLocalVector(LocalVector &lvec, const TVector &vec)
Definition: local_algebra.h:739
DebugID DID_ELEM_DISC_ASSEMBLE_UTIL("ELEM_DISC_ASSEMBLE_UTIL")
Definition: elem_disc_assemble_util.h:65
Definition: elem_disc_interface.h:77
size_t id()
Definition: util_algebra_dependent.h:130
static AlgebraTypeIDProvider & instance()
Definition: util_algebra_dependent.h:107