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  }
1669  UG_CATCH_THROW("(instationary) PrepareTimestep: Cannot create Data Evaluator.");
1670  }
1671 
1672 // //////////////////////////////////////////////////////////////////////////////
1673 // Finish Timestep (for instationary problems)
1674 // /////////////////////////////////////////////////////////////////////////////
1675 public:
1687  static void
1689  (
1690  const std::vector<IElemDisc<domain_type>*>& vElemDisc,
1692  bool bNonRegularGrid,
1695  )
1696  {
1697  // get current time and vector
1698  const vector_type& u = *vSol->solution(0);
1699 
1700  // create local time series
1701  LocalVectorTimeSeries locTimeSeries;
1702  locTimeSeries.read_times(vSol);
1703 
1704  try
1705  {
1707  vElemDisc, dd->function_pattern(), bNonRegularGrid,
1708  &locTimeSeries);
1709  Eval.set_time_point(0);
1710 
1711  // finish time step
1712  try
1713  {
1715  size_t algebra_index = bridge::AlgebraTypeIDProvider::instance().id<algebra_type>();
1716  Eval.finish_timestep(vSol->time(0), &vp, algebra_index);
1717  }
1718  UG_CATCH_THROW("(instationary) FinishTimestep: Cannot prepare time step.");
1719 
1720  }
1721  UG_CATCH_THROW("(instationary) FinishTimestep: Cannot create Data Evaluator.");
1722  }
1723 
1725 // Finish Timestep Elem (for instationary problems)
1727 
1728 public:
1744  template <typename TElem, typename TIterator>
1745  static void
1746  FinishTimestepElem(const std::vector<IElemDisc<domain_type>*>& vElemDisc,
1747  ConstSmartPtr<domain_type> spDomain,
1749  TIterator iterBegin,
1750  TIterator iterEnd,
1751  int si, bool bNonRegularGrid,
1754  {
1755  // check if there are any elements at all, otherwise return immediately
1756  if(iterBegin == iterEnd) return;
1757 
1758  // reference object id
1760 
1761  // storage for corner coordinates
1763 
1764  // get current time and vector
1765  const vector_type& u = *vSol->solution(0);
1766 
1767  // create local time series
1768  LocalVectorTimeSeries locTimeSeries;
1769  locTimeSeries.read_times(vSol);
1770 
1771  // prepare for given elem discs
1772  try
1773  {
1775  vElemDisc, dd->function_pattern(), bNonRegularGrid,
1776  &locTimeSeries);
1777  Eval.set_time_point(0);
1778 
1779  // prepare loop
1780  Eval.prepare_elem_loop(id, si);
1781 
1782  // local algebra
1783  LocalIndices ind; LocalVector locU;
1784 
1785  // Loop over all elements
1786  for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
1787  {
1788  // get Element
1789  TElem* elem = *iter;
1790 
1791  // get corner coordinates
1792  FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
1793 
1794  // check if elem is skipped from assembling
1795  if(!spAssTuner->element_used(elem)) continue;
1796 
1797  // get global indices
1798  dd->indices(elem, ind, Eval.use_hanging());
1799 
1800  // adapt local algebra
1801  locU.resize(ind);
1802 
1803  // read local values of u
1804  GetLocalVector(locU, u);
1805 
1806  // read local values of time series
1807  if(Eval.time_series_needed())
1808  locTimeSeries.read_values(vSol, ind);
1809 
1810  // finish timestep
1811  try
1812  {
1813  Eval.finish_timestep_elem(locTimeSeries.time(0), locU, elem, vCornerCoords);
1814  }
1815  UG_CATCH_THROW("(instationary) FinishTimestepElem: Cannot finish timestep.");
1816  }
1817 
1818  }
1819  UG_CATCH_THROW("(instationary) FinishTimestepElem: Cannot create Data Evaluator");
1820  }
1821 
1823 // Error estimator (stationary)
1825 
1826 public:
1841  template <typename TElem, typename TIterator>
1842  static void
1844  (
1845  const std::vector<IElemError<domain_type>*>& vElemDisc,
1846  ConstSmartPtr<domain_type> spDomain,
1848  TIterator iterBegin,
1849  TIterator iterEnd,
1850  int si,
1851  bool bNonRegularGrid,
1852  const vector_type& u
1853  )
1854  {
1855  // check if there are any elements at all, otherwise return immediately
1856  if(iterBegin == iterEnd) return;
1857 
1858  // reference object id
1860 
1861  // storage for corner coordinates
1863 
1864  // prepare for given elem discs
1865  try
1866  {
1868  vElemDisc, dd->function_pattern(), bNonRegularGrid);
1869 
1870  // prepare element loop
1871  Eval.prepare_err_est_elem_loop(id, si);
1872 
1873  // local indices and local algebra
1874  LocalIndices ind; LocalVector locU;
1875 
1876  // Loop over all elements
1877  for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
1878  {
1879  // get Element
1880  TElem* elem = *iter;
1881 
1882  // get corner coordinates
1883  FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
1884 
1885  // get global indices
1886  dd->indices(elem, ind, Eval.use_hanging());
1887 
1888  // adapt local algebra
1889  locU.resize(ind);
1890 
1891  // read local values of u
1892  GetLocalVector(locU, u);
1893 
1894  // prepare element
1895  try
1896  {
1897  Eval.prepare_err_est_elem(locU, elem, vCornerCoords, ind, false);
1898  }
1899  UG_CATCH_THROW("(stationary) AssembleRhs: Cannot prepare element.");
1900 
1901  // assemble stiffness part
1902  try
1903  {
1904  Eval.compute_err_est_A_elem(locU, elem, vCornerCoords, ind);
1905  }
1906  UG_CATCH_THROW("AssembleErrorEstimator: Cannot assemble the error estimator for stiffness part.");
1907 
1908  // assemble rhs part
1909  try
1910  {
1911  Eval.compute_err_est_rhs_elem(elem, vCornerCoords, ind);
1912  }
1913  UG_CATCH_THROW("AssembleErrorEstimator: Cannot assemble the error estimator for stiffness part.");
1914 
1915  }
1916 
1917  // finish element loop
1918  try
1919  {
1920  Eval.finish_err_est_elem_loop();
1921  }
1922  UG_CATCH_THROW("AssembleErrorEstimator: Cannot finish element loop.");
1923 
1924  }
1925  UG_CATCH_THROW("AssembleErrorEstimator: Cannot create Data Evaluator.");
1926  }
1927 
1929 // Error estimator (for time-dependent problems)
1931 
1932 public:
1949  template <typename TElem, typename TIterator>
1950  static void
1952  (
1953  const std::vector<IElemError<domain_type>*>& vElemDisc,
1954  ConstSmartPtr<domain_type> spDomain,
1956  TIterator iterBegin,
1957  TIterator iterEnd,
1958  int si,
1959  bool bNonRegularGrid,
1960  const std::vector<number>& vScaleMass,
1961  const std::vector<number>& vScaleStiff,
1963  )
1964  {
1965  // check if there are any elements at all, otherwise return immediately
1966  if (iterBegin == iterEnd) return;
1967 
1968  // reference object id
1970 
1971  // storage for corner coordinates
1973 
1974  // check time scheme
1975  if(vScaleMass.size() != vScaleStiff.size())
1976  UG_THROW("(instationary) AssembleErrorEstimator: s_a and s_m must have same size.");
1977 
1978  if(vSol->size() < vScaleStiff.size())
1979  UG_THROW("(instationary) AssembleErrorEstimator: Time stepping scheme needs at "
1980  "least "<<vScaleStiff.size()<<" time steps, but only "<<
1981  vSol->size() << " passed.");
1982 
1983  // create local time series
1984  LocalVectorTimeSeries locTimeSeries;
1985  locTimeSeries.read_times(vSol);
1986 
1987  // prepare for given elem discs
1988  try
1989  {
1991  vElemDisc, dd->function_pattern(), bNonRegularGrid,
1992  &locTimeSeries, &vScaleMass, &vScaleStiff);
1993 
1994  // prepare element loop
1995  Eval.prepare_err_est_elem_loop(id, si);
1996 
1997  // local indices and local algebra
1998  LocalIndices ind;
1999 
2000  // loop over all elements
2001  for (TIterator iter = iterBegin; iter != iterEnd; ++iter)
2002  {
2003  // get Element
2004  TElem* elem = *iter;
2005 
2006  // get corner coordinates
2007  FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
2008 
2009  // get global indices
2010  dd->indices(elem, ind, Eval.use_hanging());
2011 
2012  // read local values of time series
2013  locTimeSeries.read_values(vSol, ind);
2014 
2015  // loop all time points and assemble them
2016  for (std::size_t t = 0; t < vScaleStiff.size(); ++t)
2017  {
2018  // get local solution at timepoint
2019  LocalVector& locU = locTimeSeries.solution(t);
2020  Eval.set_time_point(t);
2021 
2022  // prepare element
2023  try
2024  {
2025  Eval.prepare_err_est_elem(locU, elem, vCornerCoords, ind, false);
2026  }
2027  UG_CATCH_THROW("AssembleErrorEstimator: Cannot prepare element.");
2028 
2029  // assemble stiffness part
2030  try
2031  {
2032  Eval.compute_err_est_A_elem(locU, elem, vCornerCoords, ind, vScaleMass[t], vScaleStiff[t]);
2033  }
2034  UG_CATCH_THROW("AssembleErrorEstimator: Cannot assemble the error estimator for stiffness part.");
2035 
2036  // assemble mass part
2037  try
2038  {
2039  Eval.compute_err_est_M_elem(locU, elem, vCornerCoords, ind, vScaleMass[t], vScaleStiff[t]);
2040  }
2041  UG_CATCH_THROW("AssembleErrorEstimator: Cannot assemble the error estimator for stiffness part.");
2042 
2043  // assemble rhs part
2044  try
2045  {
2046  Eval.compute_err_est_rhs_elem(elem, vCornerCoords, ind, vScaleMass[t], vScaleStiff[t]);
2047  }
2048  UG_CATCH_THROW("AssembleErrorEstimator: Cannot assemble the error estimator for stiffness part.");
2049  }
2050  }
2051 
2052  // finish element loop
2053  try
2054  {
2055  Eval.finish_err_est_elem_loop();
2056  }
2057  UG_CATCH_THROW("AssembleErrorEstimator: Cannot finish element loop.");
2058 
2059  }
2060  UG_CATCH_THROW("AssembleErrorEstimator: Cannot create Data Evaluator.");
2061  }
2062 
2063 }; // class StdGlobAssembler
2064 
2065 } // end namespace ug
2066 
2067 
2068 #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:1689
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:1844
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:1952
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:1746
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