ug4
dom_disc_embb_impl.h
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2020: G-CSC, Goethe University Frankfurt
3  * Author: Dmitry Logashenko
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 /*
34  * Implementation of functions from embass.h
35  */
36 
37 namespace ug {
38 
39 /*-------- class LSGFGlobAssembler --------*/
40 
57 template <typename TDomain, typename TAlgebra, typename TExtrapolation>
58 template <typename TElem, typename TIterator>
60 AssembleStiffnessMatrix( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
63  TIterator iterBegin,
64  TIterator iterEnd,
65  int si, bool bNonRegularGrid,
66  matrix_type& A,
67  const vector_type& u,
69 {
70 // check the dof distribution for the ghost-fluid method
71  m_extrapol.check_dd (dd);
72 
73 // check if there are any elements at all, otherwise return immediately
74  if(iterBegin == iterEnd) return;
75 
76 // reference object id
78 
79 // storage for corner coordinates
81 
82  try
83  {
85  vElemDisc, dd->function_pattern(), bNonRegularGrid);
86 
87 // prepare element loop
88  Eval.prepare_elem_loop(id, si);
89 
90 // local indices and local algebra
91  LocalIndices ind; LocalVector locU; LocalMatrix locA;
92 
93 // grid level of all the elements
94  int g_level = dd->grid_level().level ();
95  if (g_level == GridLevel::TOP)
96  g_level = dd->multi_grid()->top_level ();
97 
98 // Loop over all elements
99  for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
100  {
101  // get Element
102  TElem* elem = *iter;
103 
104  // check if elem is skipped from assembling
105  if(!spAssTuner->element_used(elem)) continue;
106 
107  // get corner coordinates
108  FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
109 
110  // check whether we are inside
111  int elem_status = m_extrapol.check_elem_lsf
112  (TElem::NUM_VERTICES, elem, si, g_level, Eval.use_hanging(), vCornerCoords, 0);
113  if (elem_status < 0) continue; // this is an outer element, do not assemble it at all
114  if (m_bAssembleOnlyCut && elem_status > 0) continue; // exclude elements that are not cut
115 
116  // get global indices
117  dd->indices(elem, ind, Eval.use_hanging());
118 
119  // adapt local algebra
120  locU.resize(ind); locA.resize(ind);
121 
122  // assemble depending on whether this is an inner or an interface element
123  if (elem_status > 0) // an inner element, usual assembling
124  {
125  // read local values of u
126  GetLocalVector(locU, u);
127 
128  // prepare element
129  try
130  {
131  Eval.prepare_elem(locU, elem, id, vCornerCoords, ind, true);
132  }
133  UG_CATCH_THROW("AssembleStiffnessMatrix Cannot prepare element.");
134 
135  // Assemble JA
136  locA = 0.0;
137  try
138  {
139  Eval.add_jac_A_elem(locA, locU, elem, vCornerCoords);
140  }
141  UG_CATCH_THROW("AssembleStiffnessMatrix: Cannot compute Jacobian (A).");
142 
143  // send local to global matrix
144  try{
145  spAssTuner->add_local_mat_to_global(A, locA, dd);
146  }
147  UG_CATCH_THROW("AssembleStiffnessMatrix: Cannot add local matrix.");
148  }
149  else
150  {
151  // assemble the matrix for every corner separately
152  for (size_t base_co = 0; base_co < TElem::NUM_VERTICES; base_co++)
153  if (m_extrapol.corner_inside (base_co)) /* ... only for the corners 'inside' */
154  {
155  // read local values of u
156  GetLocalVector(locU, u);
157 
158  // extrapolate the inner values to the outer ones
159  m_extrapol.template extrapolate_sol_by_lsf<TElem> (locU, base_co);
160 
161  // prepare element
162  try
163  {
164  Eval.prepare_elem(locU, elem, id, vCornerCoords, ind, true);
165  }
166  UG_CATCH_THROW("AssembleStiffnessMatrix Cannot prepare element.");
167 
168  // Assemble JA
169  locA = 0.0;
170  try
171  {
172  Eval.add_jac_A_elem(locA, locU, elem, vCornerCoords);
173  }
174  UG_CATCH_THROW("AssembleStiffnessMatrix: Cannot compute Jacobian (A).");
175 
176  // eliminate the outer connections in the local matrix
177  m_extrapol.template eliminate_extrapolated<TElem> (locA, base_co);
178 
179  // send local to global matrix
180  try{
181  spAssTuner->add_local_mat_to_global(A, locA, dd);
182  }
183  UG_CATCH_THROW("AssembleStiffnessMatrix: Cannot add local matrix.");
184  }
185  }
186  }
187 
188 // finish element loop
189  try
190  {
191  Eval.finish_elem_loop();
192  }
193  UG_CATCH_THROW("AssembleStiffnessMatrix: Cannot finish element loop.");
194 
195  }
196  UG_CATCH_THROW("AssembleStiffnessMatrix': Cannot create Data Evaluator.");
197 }
198 
215 template <typename TDomain, typename TAlgebra, typename TExtrapolation>
216 template <typename TElem, typename TIterator>
218 AssembleMassMatrix( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
221  TIterator iterBegin,
222  TIterator iterEnd,
223  int si, bool bNonRegularGrid,
224  matrix_type& M,
225  const vector_type& u,
227 {
228 // check the dof distribution for the ghost-fluid method
229  m_extrapol.check_dd (dd);
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 // grid level of all the elements
253  int g_level = dd->grid_level().level ();
254  if (g_level == GridLevel::TOP)
255  g_level = dd->multi_grid()->top_level ();
256 
257 // Loop over all elements
258  for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
259  {
260  // get Element
261  TElem* elem = *iter;
262 
263  // check if elem is skipped from assembling
264  if(!spAssTuner->element_used(elem)) continue;
265 
266  // get corner coordinates
267  FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
268 
269  // check whether we are inside
270  int elem_status = m_extrapol.check_elem_lsf
271  (TElem::NUM_VERTICES, elem, si, g_level, Eval.use_hanging(), vCornerCoords, 0);
272  if (elem_status < 0) continue; // this is an outer element, do not assemble it at all
273  if (m_bAssembleOnlyCut && elem_status > 0) continue; // exclude elements that are not cut
274 
275  // get global indices
276  dd->indices(elem, ind, Eval.use_hanging());
277 
278  // adapt local algebra
279  locU.resize(ind); locM.resize(ind);
280 
281  // assemble depending on whether this is an inner or an interface element
282  if (elem_status > 0) // an inner element, usual assembling
283  {
284  // read local values of u
285  GetLocalVector(locU, u);
286 
287  // prepare element
288  try
289  {
290  Eval.prepare_elem(locU, elem, id, vCornerCoords, ind, true);
291  }
292  UG_CATCH_THROW("AssembleMassMatrix: Cannot prepare element.");
293 
294  // Assemble JM
295  locM = 0.0;
296  try
297  {
298  Eval.add_jac_M_elem(locM, locU, elem, vCornerCoords);
299  }
300  UG_CATCH_THROW("AssembleMassMatrix: Cannot compute Jacobian (M).");
301 
302  // send local to global matrix
303  try{
304  spAssTuner->add_local_mat_to_global(M, locM, dd);
305  }
306  UG_CATCH_THROW("AssembleMassMatrix: Cannot add local matrix.");
307  }
308  else // an interface element
309  {
310  // assemble the matrix for every corner separately
311  for (size_t base_co = 0; base_co < TElem::NUM_VERTICES; base_co++)
312  if (m_extrapol.corner_inside (base_co)) /* ... only for the corners 'inside' */
313  {
314  // read local values of u
315  GetLocalVector(locU, u);
316 
317  // extrapolate the inner values to the outer ones
318  m_extrapol.template extrapolate_sol_by_lsf<TElem> (locU, base_co);
319 
320  // prepare element
321  try
322  {
323  Eval.prepare_elem(locU, elem, id, vCornerCoords, ind, true);
324  }
325  UG_CATCH_THROW("AssembleMassMatrix: Cannot prepare element.");
326 
327  // Assemble JM
328  locM = 0.0;
329  try
330  {
331  Eval.add_jac_M_elem(locM, locU, elem, vCornerCoords);
332  }
333  UG_CATCH_THROW("AssembleMassMatrix: Cannot compute Jacobian (M).");
334 
335  // eliminate the outer connections in the local matrix
336  m_extrapol.template eliminate_extrapolated<TElem> (locM, base_co);
337 
338  // send local to global matrix
339  try{
340  spAssTuner->add_local_mat_to_global(M, locM, dd);
341  }
342  UG_CATCH_THROW("AssembleMassMatrix: Cannot add local matrix.");
343  }
344  }
345  }
346 
347 // finish element loop
348  try
349  {
350  Eval.finish_elem_loop();
351  }
352  UG_CATCH_THROW("AssembleMassMatrix: Cannot finish element loop.");
353 
354  }
355  UG_CATCH_THROW("AssembleMassMatrix: Cannot create Data Evaluator.");
356 }
357 
373 template <typename TDomain, typename TAlgebra, typename TExtrapolation>
374 template <typename TElem, typename TIterator>
376 AssembleJacobian( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
379  TIterator iterBegin,
380  TIterator iterEnd,
381  int si, bool bNonRegularGrid,
382  matrix_type& J,
383  const vector_type& u,
385 {
386 // check the dof distribution for the ghost-fluid method
387  m_extrapol.check_dd (dd);
388 
389 // check if there are any elements at all, otherwise return immediately
390  if(iterBegin == iterEnd) return;
391 
392 // reference object id
394 
395 // storage for corner coordinates
397 
398 // prepare for given elem discs
399  try
400  {
402  vElemDisc, dd->function_pattern(), bNonRegularGrid);
403 
404 // prepare element loop
405  Eval.prepare_elem_loop(id, si);
406 
407 // local indices and local algebra
408  LocalIndices ind; LocalVector locU; LocalMatrix locJ;
409 
410 // grid level of all the elements
411  int g_level = dd->grid_level().level ();
412  if (g_level == GridLevel::TOP)
413  g_level = dd->multi_grid()->top_level ();
414 
415 // Loop over all elements
416  for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
417  {
418  // get Element
419  TElem* elem = *iter;
420 
421  // check if elem is skipped from assembling
422  if(!spAssTuner->element_used(elem)) continue;
423 
424  // get corner coordinates
425  FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
426 
427  // check whether we are inside
428  int elem_status = m_extrapol.check_elem_lsf
429  (TElem::NUM_VERTICES, elem, si, g_level, Eval.use_hanging(), vCornerCoords, 0);
430  if (elem_status < 0) continue; // this is an outer element, do not assemble it at all
431  if (m_bAssembleOnlyCut && elem_status > 0) continue; // exclude elements that are not cut
432 
433  // get global indices
434  dd->indices(elem, ind, Eval.use_hanging());
435 
436  // adapt local algebra
437  locU.resize(ind); locJ.resize(ind);
438 
439  // assemble depending on whether this is an inner or an interface element
440  if (elem_status > 0) // an inner element, usual assembling
441  {
442  // read local values of u
443  GetLocalVector(locU, u);
444 
445  // prepare element
446  try
447  {
448  Eval.prepare_elem(locU, elem, id, vCornerCoords, ind, true);
449  }
450  UG_CATCH_THROW("(stationary) AssembleJacobian: Cannot prepare element.");
451 
452  // reset local algebra
453  locJ = 0.0;
454 
455  // Assemble JA
456  try
457  {
458  Eval.add_jac_A_elem(locJ, locU, elem, vCornerCoords);
459  }
460  UG_CATCH_THROW("(stationary) AssembleJacobian: Cannot compute Jacobian (A).");
461 
462  // send local to global matrix
463  try{
464  spAssTuner->add_local_mat_to_global(J, locJ, dd);
465  }
466  UG_CATCH_THROW("(stationary) AssembleJacobian: Cannot add local matrix.");
467  }
468  else // an interface element
469  {
470  // assemble the jacobian for every corner separately
471  for (size_t base_co = 0; base_co < TElem::NUM_VERTICES; base_co++)
472  if (m_extrapol.corner_inside (base_co)) /* ... only for the corners 'inside' */
473  {
474  // read local values of u
475  GetLocalVector(locU, u);
476 
477  // extrapolate the inner values to the outer ones
478  m_extrapol.template extrapolate_sol_by_lsf<TElem> (locU, base_co);
479 
480  // prepare element
481  try
482  {
483  Eval.prepare_elem(locU, elem, id, vCornerCoords, ind, true);
484  }
485  UG_CATCH_THROW("(stationary) AssembleJacobian: Cannot prepare element.");
486 
487  // reset local algebra
488  locJ = 0.0;
489 
490  // Assemble JA
491  try
492  {
493  Eval.add_jac_A_elem(locJ, locU, elem, vCornerCoords);
494  }
495  UG_CATCH_THROW("(stationary) AssembleJacobian: Cannot compute Jacobian (A).");
496 
497  // eliminate the outer connections in the local matrix
498  m_extrapol.template eliminate_extrapolated<TElem> (locJ, base_co);
499 
500  // send local to global matrix
501  try{
502  spAssTuner->add_local_mat_to_global(J, locJ, dd);
503  }
504  UG_CATCH_THROW("(stationary) AssembleJacobian: Cannot add local matrix.");
505  }
506  }
507  }
508 
509 // finish element loop
510  try
511  {
512  Eval.finish_elem_loop();
513  }
514  UG_CATCH_THROW("(stationary) AssembleJacobian: Cannot finish element loop.");
515 
516  }
517  UG_CATCH_THROW("(stationary) AssembleJacobian: Cannot create Data Evaluator.");
518 }
519 
538 template <typename TDomain, typename TAlgebra, typename TExtrapolation>
539 template <typename TElem, typename TIterator>
541 AssembleJacobian( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
544  TIterator iterBegin,
545  TIterator iterEnd,
546  int si, bool bNonRegularGrid,
547  matrix_type& J,
549  number s_a0,
551 {
552 // check the dof distribution for the ghost-fluid method
553  m_extrapol.check_dd (dd);
554 
555 // check if there are any elements at all, otherwise return immediately
556  if(iterBegin == iterEnd) return;
557 
558 // reference object id
560 
561 // storage for corner coordinates
563 
564 // get current time and vector
565  const vector_type& u = *vSol->solution(0);
566 
567 // create local time series
568  LocalVectorTimeSeries locTimeSeries;
569  locTimeSeries.read_times(vSol);
570 
571 // prepare for given elem discs
572  try
573  {
575  vElemDisc, dd->function_pattern(), bNonRegularGrid,
576  &locTimeSeries);
577  Eval.set_time_point(0);
578 
579 // prepare element loop
580  Eval.prepare_elem_loop(id, si);
581 
582 // local algebra
583  LocalIndices ind; LocalVector locU; LocalMatrix locJ;
584 
585 // grid level of all the elements
586  int g_level = dd->grid_level().level ();
587  if (g_level == GridLevel::TOP)
588  g_level = dd->multi_grid()->top_level ();
589 
590 // Loop over all elements
591  for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
592  {
593  // get Element
594  TElem* elem = *iter;
595 
596  // check if elem is skipped from assembling
597  if(!spAssTuner->element_used(elem)) continue;
598 
599  // get corner coordinates
600  FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
601 
602  // check whether we are inside
603  int elem_status = m_extrapol.check_elem_lsf
604  (TElem::NUM_VERTICES, elem, si, g_level, Eval.use_hanging(), vCornerCoords, vSol->time(0));
605  if (elem_status < 0) continue; // this is an outer element, do not assemble it at all
606  if (m_bAssembleOnlyCut && elem_status > 0) continue; // exclude elements that are not cut
607 
608  // get global indices
609  dd->indices(elem, ind, Eval.use_hanging());
610 
611  // adapt local algebra
612  locU.resize(ind); locJ.resize(ind);
613 
614  // assemble depending on whether this is an inner or an interface element
615  if (elem_status > 0) // an inner element, usual assembling
616  {
617  // read local values of u
618  GetLocalVector(locU, u);
619 
620  // read local values of time series
621  if(Eval.time_series_needed())
622  locTimeSeries.read_values(vSol, ind);
623 
624  // prepare element
625  try
626  {
627  Eval.prepare_elem(locU, elem, id, vCornerCoords, ind, true);
628  }
629  UG_CATCH_THROW("(instationary) AssembleJacobian: Cannot prepare element.");
630 
631  // reset local algebra
632  locJ = 0.0;
633 
634  // Assemble JA
635  try
636  {
637  Eval.add_jac_A_elem(locJ, locU, elem, vCornerCoords, PT_INSTATIONARY);
638  locJ *= s_a0;
639 
640  Eval.add_jac_A_elem(locJ, locU, elem, vCornerCoords, PT_STATIONARY);
641  }
642  UG_CATCH_THROW("(instationary) AssembleJacobian: Cannot compute Jacobian (A).");
643 
644  // Assemble JM
645  try
646  {
647  Eval.add_jac_M_elem(locJ, locU, elem, vCornerCoords, PT_INSTATIONARY);
648  }
649  UG_CATCH_THROW("(instationary) AssembleJacobian: Cannot compute Jacobian (M).");
650 
651  // send local to global matrix
652  try{
653  spAssTuner->add_local_mat_to_global(J, locJ, dd);
654  }
655  UG_CATCH_THROW("(instationary) AssembleJacobian: Cannot add local matrix.");
656  }
657  else // an interface element
658  {
659  // assemble the jacobian for every corner separately
660  for (size_t base_co = 0; base_co < TElem::NUM_VERTICES; base_co++)
661  if (m_extrapol.corner_inside (base_co)) /* ... only for the corners 'inside' */
662  {
663  // read local values of u
664  GetLocalVector(locU, u);
665 
666  // extrapolate the inner values to the outer ones
667  m_extrapol.template extrapolate_sol_by_lsf<TElem> (locU, base_co);
668 
669  // read local values of time series
670  if(Eval.time_series_needed())
671  {
672  locTimeSeries.read_values(vSol, ind);
673  for (size_t i = 0; i < locTimeSeries.size (); i++)
674  m_extrapol.template extrapolate_sol_by_lsf<TElem>
675  (locTimeSeries.solution (i), base_co);
676  }
677 
678  // prepare element
679  try
680  {
681  Eval.prepare_elem(locU, elem, id, vCornerCoords, ind, true);
682  }
683  UG_CATCH_THROW("(instationary) AssembleJacobian: Cannot prepare element.");
684 
685  // reset local algebra
686  locJ = 0.0;
687 
688  // Assemble JA
689  try
690  {
691  Eval.add_jac_A_elem(locJ, locU, elem, vCornerCoords, PT_INSTATIONARY);
692  locJ *= s_a0;
693 
694  Eval.add_jac_A_elem(locJ, locU, elem, vCornerCoords, PT_STATIONARY);
695  }
696  UG_CATCH_THROW("(instationary) AssembleJacobian: Cannot compute Jacobian (A).");
697 
698  // Assemble JM
699  try
700  {
701  Eval.add_jac_M_elem(locJ, locU, elem, vCornerCoords, PT_INSTATIONARY);
702  }
703  UG_CATCH_THROW("(instationary) AssembleJacobian: Cannot compute Jacobian (M).");
704 
705  // eliminate the outer connections in the local matrix
706  m_extrapol.template eliminate_extrapolated<TElem> (locJ, base_co);
707 
708  // send local to global matrix
709  try{
710  spAssTuner->add_local_mat_to_global(J, locJ, dd);
711  }
712  UG_CATCH_THROW("(instationary) AssembleJacobian: Cannot add local matrix.");
713  }
714  }
715  }
716 
717 // finish element loop
718  try
719  {
720  Eval.finish_elem_loop();
721  }
722  UG_CATCH_THROW("(instationary) AssembleJacobian: Cannot finish element loop.");
723 
724  }
725  UG_CATCH_THROW("(instationary) AssembleJacobian: Cannot create Data Evaluator.");
726 }
727 
744 template <typename TDomain, typename TAlgebra, typename TExtrapolation>
745 template <typename TElem, typename TIterator>
747 AssembleDefect( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
750  TIterator iterBegin,
751  TIterator iterEnd,
752  int si, bool bNonRegularGrid,
753  vector_type& d,
754  const vector_type& u,
756 {
757 // check the dof distribution for the ghost-fluid method
758  m_extrapol.check_dd (dd);
759 
760 // check if at least one element exists, else return
761  if(iterBegin == iterEnd) return;
762 
763 // reference object id
765 
766 // storage for corner coordinates
768 
769 // prepare for given elem discs
770  try
771  {
773  vElemDisc, dd->function_pattern(), bNonRegularGrid);
774 
775 // prepare element loop
776  Eval.prepare_elem_loop(id, si);
777 
778 // local indices and local algebra
779  LocalIndices ind; LocalVector locU, locD, tmpLocD;
780 
781 // grid level of all the elements
782  int g_level = dd->grid_level().level ();
783  if (g_level == GridLevel::TOP)
784  g_level = dd->multi_grid()->top_level ();
785 
786 // Loop over all elements
787  for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
788  {
789  // get Element
790  TElem* elem = *iter;
791 
792  // check if elem is skipped from assembling
793  if(!spAssTuner->element_used(elem)) continue;
794 
795  // get corner coordinates
796  FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
797 
798  // check whether we are inside
799  int elem_status = m_extrapol.check_elem_lsf
800  (TElem::NUM_VERTICES, elem, si, g_level, Eval.use_hanging(), vCornerCoords, 0);
801  if (elem_status < 0) continue; // this is an outer element, do not assemble it at all
802  if (m_bAssembleOnlyCut && elem_status > 0) continue; // exclude elements that are not cut
803 
804  // get global indices
805  dd->indices(elem, ind, Eval.use_hanging());
806 
807  // adapt local algebra
808  locU.resize(ind); locD.resize(ind); tmpLocD.resize(ind);
809 
810  // assemble depending on whether this is an inner or an interface element
811  if (elem_status > 0) // an inner element, usual assembling
812  {
813  // read local values of u
814  GetLocalVector(locU, u);
815 
816  // prepare element
817  try
818  {
819  Eval.prepare_elem(locU, elem, id, vCornerCoords, ind);
820  }
821  UG_CATCH_THROW("(stationary) AssembleDefect: Cannot prepare element.");
822 
823  // ANALOG to 'domain_disc_elem()' - modifies the solution, used
824  // for computing the defect
825  if( spAssTuner->modify_solution_enabled() )
826  {
827  LocalVector& modLocU = locU;
828  try{
829  spAssTuner->modify_LocalSol(modLocU, locU, dd);
830  } UG_CATCH_THROW("Cannot modify local solution.");
831 
832  // recopy modified LocalVector:
833  locU = modLocU;
834  }
835 
836  // reset local algebra
837  locD = 0.0;
838 
839  // Assemble A
840  try
841  {
842  Eval.add_def_A_elem(locD, locU, elem, vCornerCoords);
843  }
844  UG_CATCH_THROW("(stationary) AssembleDefect: Cannot compute Defect (A).");
845 
846  // Assemble rhs
847  try
848  {
849  tmpLocD = 0.0;
850  Eval.add_rhs_elem(tmpLocD, elem, vCornerCoords);
851  locD.scale_append(-1, tmpLocD);
852 
853  }
854  UG_CATCH_THROW("(stationary) AssembleDefect: Cannot compute Rhs.");
855 
856  // send local to global defect
857  try{
858  spAssTuner->add_local_vec_to_global(d, locD, dd);
859  }
860  UG_CATCH_THROW("(stationary) AssembleDefect: Cannot add local vector.");
861  }
862  else // an interface element
863  {
864  // assemble the defect for every corner separately
865  for (size_t base_co = 0; base_co < TElem::NUM_VERTICES; base_co++)
866  if (m_extrapol.corner_inside (base_co)) /* ... only for the corners 'inside' */
867  {
868  // read local values of u
869  GetLocalVector(locU, u);
870 
871  // extrapolate the inner values to the outer ones
872  m_extrapol.template extrapolate_sol_by_lsf<TElem> (locU, base_co);
873 
874  // prepare element
875  try
876  {
877  Eval.prepare_elem(locU, elem, id, vCornerCoords, ind);
878  }
879  UG_CATCH_THROW("(stationary) AssembleDefect: Cannot prepare element.");
880 
881  // ANALOG to 'domain_disc_elem()' - modifies the solution, used
882  // for computing the defect
883  if( spAssTuner->modify_solution_enabled() )
884  {
885  LocalVector& modLocU = locU;
886  try{
887  spAssTuner->modify_LocalSol(modLocU, locU, dd);
888  } UG_CATCH_THROW("Cannot modify local solution.");
889 
890  // recopy modified LocalVector:
891  locU = modLocU;
892  }
893 
894  // reset local algebra
895  locD = 0.0;
896 
897  // Assemble A
898  try
899  {
900  Eval.add_def_A_elem(locD, locU, elem, vCornerCoords);
901  }
902  UG_CATCH_THROW("(stationary) AssembleDefect: Cannot compute Defect (A).");
903 
904  // Assemble rhs
905  try
906  {
907  tmpLocD = 0.0;
908  Eval.add_rhs_elem(tmpLocD, elem, vCornerCoords);
909  locD.scale_append(-1, tmpLocD);
910 
911  }
912  UG_CATCH_THROW("(stationary) AssembleDefect: Cannot compute Rhs.");
913 
914  // eliminate the outer connections in the local matrix
915  m_extrapol.template clear_outer_vectors<TElem> (locD, base_co);
916 
917  // send local to global defect
918  try{
919  spAssTuner->add_local_vec_to_global(d, locD, dd);
920  }
921  UG_CATCH_THROW("(stationary) AssembleDefect: Cannot add local vector.");
922  }
923  }
924  }
925 
926 // finish element loop
927  try
928  {
929  Eval.finish_elem_loop();
930  }
931  UG_CATCH_THROW("(stationary) AssembleDefect: Cannot finish element loop.");
932 
933  }
934  UG_CATCH_THROW("(stationary) AssembleDefect: Cannot create Data Evaluator.");
935 }
936 
955 template <typename TDomain, typename TAlgebra, typename TExtrapolation>
956 template <typename TElem, typename TIterator>
958 AssembleDefect( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
961  TIterator iterBegin,
962  TIterator iterEnd,
963  int si, bool bNonRegularGrid,
964  vector_type& d,
966  const std::vector<number>& vScaleMass,
967  const std::vector<number>& vScaleStiff,
969 {
970 // check the dof distribution for the ghost-fluid method
971  m_extrapol.check_dd (dd);
972 
973 // check if there are any elements at all, otherwise return immediately
974  if(iterBegin == iterEnd) return;
975 
976 // reference object id
978 
979 // storage for corner coordinates
981 
982 // check time scheme
983  if(vScaleMass.size() != vScaleStiff.size())
984  UG_THROW("(instationary) AssembleDefect: s_a and s_m must have same size.");
985 
986  if(vSol->size() < vScaleStiff.size())
987  UG_THROW("(instationary) AssembleDefect: Time stepping scheme needs at "
988  "least "<<vScaleStiff.size()<<" time steps, but only "<<
989  vSol->size() << " passed.");
990 
991 // create local time series
992  LocalVectorTimeSeries locTimeSeries;
993  locTimeSeries.read_times(vSol);
994 
995 // prepare for given elem discs
996  try
997  {
999  vElemDisc, dd->function_pattern(), bNonRegularGrid,
1000  &locTimeSeries, &vScaleMass, &vScaleStiff);
1001 
1002 // prepare element loop
1003  Eval.prepare_elem_loop(id, si);
1004 
1005 // local indices and local algebra
1006  LocalIndices ind; LocalVector locD, tmpLocD;
1007 
1008 // grid level of all the elements
1009  int g_level = dd->grid_level().level ();
1010  if (g_level == GridLevel::TOP)
1011  g_level = dd->multi_grid()->top_level ();
1012 
1013 // Loop over all elements
1014  for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
1015  {
1016  // get Element
1017  TElem* elem = *iter;
1018 
1019  // check if elem is skipped from assembling
1020  if(!spAssTuner->element_used(elem)) continue;
1021 
1022  // get corner coordinates
1023  FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
1024 
1025  // check whether we are inside
1026  int elem_status = m_extrapol.check_elem_lsf
1027  (TElem::NUM_VERTICES, elem, si, g_level, Eval.use_hanging(), vCornerCoords, vSol->time(0));
1028  if (elem_status < 0) continue; // this is an outer element, do not assemble it at all
1029  if (m_bAssembleOnlyCut && elem_status > 0) continue; // exclude elements that are not cut
1030 
1031  // get global indices
1032  dd->indices(elem, ind, Eval.use_hanging());
1033 
1034  // adapt local algebra
1035  locD.resize(ind); tmpLocD.resize(ind);
1036 
1037  // assemble depending on whether this is an inner or an interface element
1038  if (elem_status > 0) // an inner element, usual assembling
1039  {
1040  // read local values of time series
1041  locTimeSeries.read_values(vSol, ind);
1042 
1043  // reset contribution of this element
1044  locD = 0.0;
1045 
1046  // loop all time points and assemble them
1047  for(size_t t = 0; t < vScaleStiff.size(); ++t)
1048  {
1049  number scale_stiff = vScaleStiff[t];
1050 
1051  // get local solution at timepoint
1052  LocalVector& locU = locTimeSeries.solution(t);
1053  Eval.set_time_point(t);
1054 
1055  // prepare element
1056  try
1057  {
1058  Eval.prepare_elem(locU, elem, id, vCornerCoords, ind, false);
1059  }
1060  UG_CATCH_THROW("(instationary) AssembleDefect: Cannot prepare element.");
1061 
1062  // Assemble M
1063  try
1064  {
1065  tmpLocD = 0.0;
1066  Eval.add_def_M_elem(tmpLocD, locU, elem, vCornerCoords, PT_INSTATIONARY);
1067  locD.scale_append(vScaleMass[t], tmpLocD);
1068  }
1069  UG_CATCH_THROW("(instationary) AssembleDefect: Cannot compute Defect (M).");
1070 
1071  // Assemble A
1072  try
1073  {
1074  if(scale_stiff != 0.0)
1075  {
1076  tmpLocD = 0.0;
1077  Eval.add_def_A_elem(tmpLocD, locU, elem, vCornerCoords, PT_INSTATIONARY);
1078  locD.scale_append(scale_stiff, tmpLocD);
1079  }
1080 
1081  if(t == 0)
1082  Eval.add_def_A_elem(locD, locU, elem, vCornerCoords, PT_STATIONARY);
1083  }
1084  UG_CATCH_THROW("(instationary) AssembleDefect: Cannot compute Defect (A).");
1085 
1086 
1087  // Assemble defect for explicit reaction_rate, reaction and source
1088  if( t == 1 ) // only valid at lowest timediscretization order
1089  {
1090  tmpLocD = 0.0;
1091  try
1092  {
1093  Eval.add_def_A_expl_elem(tmpLocD, locU, elem, vCornerCoords, PT_INSTATIONARY);
1094  }
1095  UG_CATCH_THROW("(instationary) AssembleDefect explizit: Cannot compute Defect (A).");
1096 
1097  // UG_ASSERT(vScaleStiff.size() == 2, "Only one step method supported.");
1098  const number dt = vSol->time(0)-vSol->time(1);
1099  locD.scale_append(dt, tmpLocD);
1100  }
1101 
1102 
1103  // Assemble rhs
1104  try
1105  {
1106  if(scale_stiff != 0.0)
1107  {
1108  tmpLocD = 0.0;
1109  Eval.add_rhs_elem(tmpLocD, elem, vCornerCoords, PT_INSTATIONARY);
1110  locD.scale_append( - scale_stiff, tmpLocD);
1111  }
1112 
1113  if(t == 0)
1114  {
1115  tmpLocD = 0.0;
1116  Eval.add_rhs_elem(tmpLocD, elem, vCornerCoords, PT_STATIONARY);
1117  locD.scale_append( -1.0, tmpLocD);
1118  }
1119  }
1120  UG_CATCH_THROW("(instationary) AssembleDefect: Cannot compute Rhs.");
1121  }
1122 
1123  // send local to global defect
1124  try{
1125  spAssTuner->add_local_vec_to_global(d, locD, dd);
1126  }
1127  UG_CATCH_THROW("(instationary) AssembleDefect: Cannot add local vector.");
1128  }
1129  else // an interface element
1130  {
1131  // assemble the defect for every corner separately
1132  for (size_t base_co = 0; base_co < TElem::NUM_VERTICES; base_co++)
1133  if (m_extrapol.corner_inside (base_co)) /* ... only for the corners 'inside' */
1134  {
1135  // read local values of time series
1136  locTimeSeries.read_values(vSol, ind);
1137 
1138  // reset contribution of this element
1139  locD = 0.0;
1140 
1141  // loop all time points and assemble them
1142  for(size_t t = 0; t < vScaleStiff.size(); ++t)
1143  {
1144  number scale_stiff = vScaleStiff[t];
1145 
1146  // get local solution at timepoint
1147  LocalVector& locU = locTimeSeries.solution(t);
1148  Eval.set_time_point(t);
1149 
1150  // extrapolate the inner values to the outer ones
1151  m_extrapol.template extrapolate_sol_by_lsf<TElem> (locU, base_co);
1152 
1153  // prepare element
1154  try
1155  {
1156  Eval.prepare_elem(locU, elem, id, vCornerCoords, ind, false);
1157  }
1158  UG_CATCH_THROW("(instationary) AssembleDefect: Cannot prepare element.");
1159 
1160  // Assemble M
1161  try
1162  {
1163  tmpLocD = 0.0;
1164  Eval.add_def_M_elem(tmpLocD, locU, elem, vCornerCoords, PT_INSTATIONARY);
1165  locD.scale_append(vScaleMass[t], tmpLocD);
1166  }
1167  UG_CATCH_THROW("(instationary) AssembleDefect: Cannot compute Defect (M).");
1168 
1169  // Assemble A
1170  try
1171  {
1172  if(scale_stiff != 0.0)
1173  {
1174  tmpLocD = 0.0;
1175  Eval.add_def_A_elem(tmpLocD, locU, elem, vCornerCoords, PT_INSTATIONARY);
1176  locD.scale_append(scale_stiff, tmpLocD);
1177  }
1178 
1179  if(t == 0)
1180  Eval.add_def_A_elem(locD, locU, elem, vCornerCoords, PT_STATIONARY);
1181  }
1182  UG_CATCH_THROW("(instationary) AssembleDefect: Cannot compute Defect (A).");
1183 
1184 
1185  // Assemble defect for explicit reaction_rate, reaction and source
1186  if( t == 1 ) // only valid at lowest timediscretization order
1187  {
1188  tmpLocD = 0.0;
1189  try
1190  {
1191  Eval.add_def_A_expl_elem(tmpLocD, locU, elem, vCornerCoords, PT_INSTATIONARY);
1192  }
1193  UG_CATCH_THROW("(instationary) AssembleDefect explizit: Cannot compute Defect (A).");
1194 
1195 // UG_ASSERT(vScaleStiff.size() == 2, "Only one step method supported.");
1196  const number dt = vSol->time(0)-vSol->time(1);
1197  locD.scale_append(dt, tmpLocD);
1198  }
1199 
1200 
1201  // Assemble rhs
1202  try
1203  {
1204  if(scale_stiff != 0.0)
1205  {
1206  tmpLocD = 0.0;
1207  Eval.add_rhs_elem(tmpLocD, elem, vCornerCoords, PT_INSTATIONARY);
1208  locD.scale_append( - scale_stiff, tmpLocD);
1209  }
1210 
1211  if(t == 0)
1212  {
1213  tmpLocD = 0.0;
1214  Eval.add_rhs_elem(tmpLocD, elem, vCornerCoords, PT_STATIONARY);
1215  locD.scale_append( -1.0, tmpLocD);
1216  }
1217  }
1218  UG_CATCH_THROW("(instationary) AssembleDefect: Cannot compute Rhs.");
1219  }
1220 
1221  // eliminate the outer connections in the local matrix
1222  m_extrapol.template clear_outer_vectors<TElem>(locD, base_co);
1223 
1224  // send local to global defect
1225  try{
1226  spAssTuner->add_local_vec_to_global(d, locD, dd);
1227  }
1228  UG_CATCH_THROW("(instationary) AssembleDefect: Cannot add local vector.");
1229  }
1230  }
1231  }
1232 
1233 // finish element loop
1234  try
1235  {
1236  Eval.finish_elem_loop();
1237  }
1238  UG_CATCH_THROW("(instationary) AssembleDefect: Cannot finish element loop.");
1239 
1240  }
1241  UG_CATCH_THROW("(instationary) AssembleDefect: Cannot create Data Evaluator.");
1242 }
1243 
1261 template <typename TDomain, typename TAlgebra, typename TExtrapolation>
1262 template <typename TElem, typename TIterator>
1264 AssembleLinear( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
1265  ConstSmartPtr<domain_type> spDomain,
1267  TIterator iterBegin,
1268  TIterator iterEnd,
1269  int si, bool bNonRegularGrid,
1270  matrix_type& A,
1271  vector_type& rhs,
1273 {
1274 // check the dof distribution for the ghost-fluid method
1275  m_extrapol.check_dd (dd);
1276 
1277 // check if there are any elements at all, otherwise return immediately
1278  if(iterBegin == iterEnd) return;
1279 
1280 // reference object id
1282 
1283 // storage for corner coordinates
1285 
1286 // prepare for given elem discs
1287  try
1288  {
1290  vElemDisc, dd->function_pattern(), bNonRegularGrid);
1291 
1292 // prepare loop
1293  Eval.prepare_elem_loop(id, si);
1294 
1295 // local indices and local algebra
1296  LocalIndices ind; LocalVector locRhs; LocalMatrix locA;
1297 
1298 // grid level of all the elements
1299  int g_level = dd->grid_level().level ();
1300  if (g_level == GridLevel::TOP)
1301  g_level = dd->multi_grid()->top_level ();
1302 
1303 // Loop over all elements
1304  for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
1305  {
1306  // get Element
1307  TElem* elem = *iter;
1308 
1309  // check if elem is skipped from assembling
1310  if(!spAssTuner->element_used(elem)) continue;
1311 
1312  // get corner coordinates
1313  FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
1314 
1315  // check whether we are inside
1316  int elem_status = m_extrapol.check_elem_lsf
1317  (TElem::NUM_VERTICES, elem, si, g_level, Eval.use_hanging(), vCornerCoords, 0);
1318  if (elem_status < 0) continue; // this is an outer element, do not assemble it at all
1319  if (m_bAssembleOnlyCut && elem_status > 0) continue; // exclude elements that are not cut
1320 
1321  // get global indices
1322  dd->indices(elem, ind, Eval.use_hanging());
1323 
1324  // adapt local algebra
1325  locRhs.resize(ind); locA.resize(ind);
1326 
1327  // prepare element
1328  try
1329  {
1330  Eval.prepare_elem(locRhs, elem, id, vCornerCoords, ind, true);
1331  }
1332  UG_CATCH_THROW("(stationary) AssembleLinear: Cannot prepare element.");
1333 
1334  // reset local algebra
1335  locA = 0.0;
1336  locRhs = 0.0;
1337 
1338  // Assemble JA
1339  try
1340  {
1341  Eval.add_jac_A_elem(locA, locRhs, elem, vCornerCoords);
1342  }
1343  UG_CATCH_THROW("(stationary) AssembleLinear: Cannot compute Jacobian (A).");
1344 
1345  // Assemble rhs
1346  try
1347  {
1348  Eval.add_rhs_elem(locRhs, elem, vCornerCoords);
1349  }
1350  UG_CATCH_THROW("(stationary) AssembleLinear: Cannot compute Rhs.");
1351 
1352  // For interface elements, correct the matrix and the right-hand side
1353  if (elem_status == 0)
1354  m_extrapol.template eliminate_extrapolated<TElem> (locA, locRhs);
1355 
1356  // send local to global matrix & rhs
1357  try{
1358  spAssTuner->add_local_mat_to_global(A, locA, dd);
1359  spAssTuner->add_local_vec_to_global(rhs, locRhs, dd);
1360  }
1361  UG_CATCH_THROW("(stationary) AssembleLinear: Cannot add local vector/matrix.");
1362  }
1363 
1364 // finish element loop
1365  try
1366  {
1367  Eval.finish_elem_loop();
1368  }
1369  UG_CATCH_THROW("(stationary) AssembleLinear: Cannot finish element loop.");
1370 
1371  }
1372  UG_CATCH_THROW("(stationary) AssembleLinear: Cannot create Data Evaluator.");
1373 }
1374 
1395 template <typename TDomain, typename TAlgebra, typename TExtrapolation>
1396 template <typename TElem, typename TIterator>
1398 AssembleLinear( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
1399  ConstSmartPtr<domain_type> spDomain,
1401  TIterator iterBegin,
1402  TIterator iterEnd,
1403  int si, bool bNonRegularGrid,
1404  matrix_type& A,
1405  vector_type& rhs,
1407  const std::vector<number>& vScaleMass,
1408  const std::vector<number>& vScaleStiff,
1410 {
1411 // check the dof distribution for the ghost-fluid method
1412  m_extrapol.check_dd (dd);
1413 
1414 // check if there are any elements at all, otherwise return immediately
1415  if(iterBegin == iterEnd) return;
1416 
1417 // reference object id
1419 
1420 // storage for corner coordinates
1422 
1423 // check time scheme
1424  if(vScaleMass.size() != vScaleStiff.size())
1425  UG_THROW("(instationary) AssembleLinear: s_a and s_m must have same size.");
1426 
1427  if(vSol->size() < vScaleStiff.size())
1428  UG_THROW("(instationary) AssembleLinear: Time stepping scheme needs at "
1429  "least "<<vScaleStiff.size()<<" time steps, but only "<<
1430  vSol->size() << " passed.");
1431 
1432 // create local time solution
1433  LocalVectorTimeSeries locTimeSeries;
1434  locTimeSeries.read_times(vSol);
1435 
1436 // prepare for given elem discs
1437  try
1438  {
1440  vElemDisc, dd->function_pattern(), bNonRegularGrid,
1441  &locTimeSeries, &vScaleMass, &vScaleStiff);
1442 
1443 // prepare loop
1444  Eval.prepare_elem_loop(id, si);
1445 
1446 // local algebra
1447  LocalIndices ind; LocalVector locRhs, tmpLocRhs; LocalMatrix locA, tmpLocA;
1448 
1449 // grid level of all the elements
1450  int g_level = dd->grid_level().level ();
1451  if (g_level == GridLevel::TOP)
1452  g_level = dd->multi_grid()->top_level ();
1453 
1454 // Loop over all elements
1455  for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
1456  {
1457  // get Element
1458  TElem* elem = *iter;
1459 
1460  // check if elem is skipped from assembling
1461  if(!spAssTuner->element_used(elem)) continue;
1462 
1463  // get corner coordinates
1464  FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
1465 
1466  // check whether we are inside
1467  int elem_status = m_extrapol.check_elem_lsf
1468  (TElem::NUM_VERTICES, elem, si, g_level, Eval.use_hanging(), vCornerCoords, vSol->time(0));
1469  if (elem_status < 0) continue; // this is an outer element, do not assemble it at all
1470  if (m_bAssembleOnlyCut && elem_status > 0) continue; // exclude elements that are not cut
1471 
1472  // get global indices
1473  dd->indices(elem, ind, Eval.use_hanging());
1474 
1475  // adapt local algebra
1476  locRhs.resize(ind); tmpLocRhs.resize(ind);
1477  locA.resize(ind); tmpLocA.resize(ind);
1478 
1479  // assemble depending on whether this is an inner or an interface element
1480  if (elem_status > 0) // an inner element, usual assembling
1481  {
1482  // read local values of time series
1483  locTimeSeries.read_values(vSol, ind);
1484  Eval.set_time_point(0);
1485 
1486  // reset element contribution
1487  locA = 0.0; locRhs = 0.0;
1488 
1490  // current time step
1491 
1492  // get local solution at time point
1493  LocalVector& locU = locTimeSeries.solution(0);
1494 
1495  // prepare element
1496  try
1497  {
1498  Eval.prepare_elem(locU, elem, id, vCornerCoords, ind, true);
1499  }
1500  UG_CATCH_THROW("(instationary) AssembleLinear: Cannot prepare element.");
1501 
1502  // Assemble JM
1503  try
1504  {
1505  tmpLocA = 0.0;
1506  Eval.add_jac_M_elem(tmpLocA, locU, elem, vCornerCoords, PT_INSTATIONARY);
1507  locA.scale_append(vScaleMass[0], tmpLocA);
1508  }
1509  UG_CATCH_THROW("(instationary) AssembleLinear: Cannot compute Jacobian (M).");
1510 
1511  // Assemble JA
1512  try
1513  {
1514  tmpLocA = 0.0;
1515  Eval.add_jac_A_elem(tmpLocA, locU, elem, vCornerCoords, PT_INSTATIONARY);
1516  locA.scale_append(vScaleStiff[0], tmpLocA);
1517 
1518  Eval.add_jac_A_elem(locA, locU, elem, vCornerCoords, PT_STATIONARY);
1519  }
1520  UG_CATCH_THROW("(instationary) AssembleLinear: Cannot compute Jacobian (A).");
1521 
1522  // Assemble rhs
1523  try
1524  {
1525  tmpLocRhs = 0.0;
1526  Eval.add_rhs_elem(tmpLocRhs, elem, vCornerCoords, PT_INSTATIONARY);
1527  locRhs.scale_append(vScaleStiff[0], tmpLocRhs);
1528 
1529  Eval.add_rhs_elem(locRhs, elem, vCornerCoords, PT_STATIONARY);
1530  }
1531  UG_CATCH_THROW("(instationary) AssembleLinear: Cannot compute Rhs.");
1532 
1534  // old time steps
1535 
1536  // loop all old time points
1537  for(size_t t = 1; t < vScaleStiff.size(); ++t)
1538  {
1539  // get local solution at time point
1540  LocalVector& locU = locTimeSeries.solution(t);
1541  Eval.set_time_point(t);
1542 
1543  // prepare element
1544  try
1545  {
1546  Eval.prepare_elem(locU, elem, id, vCornerCoords, ind, false);
1547  }
1548  UG_CATCH_THROW("(instationary) AssembleLinear: Cannot prepare element.");
1549 
1550  // Assemble dM
1551  try
1552  {
1553  tmpLocRhs = 0.0;
1554  Eval.add_def_M_elem(tmpLocRhs, locU, elem, vCornerCoords, PT_INSTATIONARY);
1555  locRhs.scale_append(-vScaleMass[t], tmpLocRhs);
1556  }
1557  UG_CATCH_THROW("(instationary) AssembleLinear: Cannot compute Jacobian (M).");
1558 
1559  // Assemble dA
1560  try
1561  {
1562  tmpLocRhs = 0.0;
1563  Eval.add_def_A_elem(tmpLocRhs, locU, elem, vCornerCoords, PT_INSTATIONARY);
1564  locRhs.scale_append(-vScaleStiff[t], tmpLocRhs);
1565  }
1566  UG_CATCH_THROW("(instationary) AssembleLinear: Cannot compute Jacobian (A).");
1567 
1568  // Assemble rhs
1569  try
1570  {
1571  tmpLocRhs = 0.0;
1572  Eval.add_rhs_elem(tmpLocRhs, elem, vCornerCoords, PT_INSTATIONARY);
1573  locRhs.scale_append(vScaleStiff[t], tmpLocRhs);
1574  }
1575  UG_CATCH_THROW("(instationary) AssembleLinear: Cannot compute Rhs.");
1576  }
1577 
1578  // send local to global matrix & rhs
1579  try{
1580  spAssTuner->add_local_mat_to_global(A, locA, dd);
1581  spAssTuner->add_local_vec_to_global(rhs, locRhs, dd);
1582  }
1583  UG_CATCH_THROW("(instationary) AssembleLinear: Cannot add local vector/matrix.");
1584  }
1585  else // an interface element
1586  {
1587  // assemble the system for every corner separately
1588  for (size_t base_co = 0; base_co < TElem::NUM_VERTICES; base_co++)
1589  if (m_extrapol.corner_inside (base_co)) /* ... only for the corners 'inside' */
1590  {
1591  // read local values of time series
1592  locTimeSeries.read_values(vSol, ind);
1593  Eval.set_time_point(0);
1594 
1595  // reset element contribution
1596  locA = 0.0; locRhs = 0.0;
1597 
1599  // current time step
1600 
1601  // get local solution at time point
1602  LocalVector& locU = locTimeSeries.solution(0);
1603 
1604  // extrapolate the inner values to the outer ones
1605  m_extrapol.template extrapolate_sol_by_lsf<TElem> (locU, base_co);
1606 
1607  // prepare element
1608  try
1609  {
1610  Eval.prepare_elem(locU, elem, id, vCornerCoords, ind, true);
1611  }
1612  UG_CATCH_THROW("(instationary) AssembleLinear: Cannot prepare element.");
1613 
1614  // Assemble JM
1615  try
1616  {
1617  tmpLocA = 0.0;
1618  Eval.add_jac_M_elem(tmpLocA, locU, elem, vCornerCoords, PT_INSTATIONARY);
1619  locA.scale_append(vScaleMass[0], tmpLocA);
1620  }
1621  UG_CATCH_THROW("(instationary) AssembleLinear: Cannot compute Jacobian (M).");
1622 
1623  // Assemble JA
1624  try
1625  {
1626  tmpLocA = 0.0;
1627  Eval.add_jac_A_elem(tmpLocA, locU, elem, vCornerCoords, PT_INSTATIONARY);
1628  locA.scale_append(vScaleStiff[0], tmpLocA);
1629 
1630  Eval.add_jac_A_elem(locA, locU, elem, vCornerCoords, PT_STATIONARY);
1631  }
1632  UG_CATCH_THROW("(instationary) AssembleLinear: Cannot compute Jacobian (A).");
1633 
1634  // Assemble rhs
1635  try
1636  {
1637  tmpLocRhs = 0.0;
1638  Eval.add_rhs_elem(tmpLocRhs, elem, vCornerCoords, PT_INSTATIONARY);
1639  locRhs.scale_append(vScaleStiff[0], tmpLocRhs);
1640 
1641  Eval.add_rhs_elem(locRhs, elem, vCornerCoords, PT_STATIONARY);
1642  }
1643  UG_CATCH_THROW("(instationary) AssembleLinear: Cannot compute Rhs.");
1644 
1646  // old time steps
1647 
1648  // loop all old time points
1649  for(size_t t = 1; t < vScaleStiff.size(); ++t)
1650  {
1651  // get local solution at time point
1652  LocalVector& locU = locTimeSeries.solution(t);
1653  Eval.set_time_point(t);
1654 
1655  // extrapolate the inner values to the outer ones
1656  m_extrapol.template extrapolate_sol_by_lsf<TElem> (locU, base_co);
1657 
1658  // prepare element
1659  try
1660  {
1661  Eval.prepare_elem(locU, elem, id, vCornerCoords, ind, false);
1662  }
1663  UG_CATCH_THROW("(instationary) AssembleLinear: Cannot prepare element.");
1664 
1665  // Assemble dM
1666  try
1667  {
1668  tmpLocRhs = 0.0;
1669  Eval.add_def_M_elem(tmpLocRhs, locU, elem, vCornerCoords, PT_INSTATIONARY);
1670  locRhs.scale_append(-vScaleMass[t], tmpLocRhs);
1671  }
1672  UG_CATCH_THROW("(instationary) AssembleLinear: Cannot compute Jacobian (M).");
1673 
1674  // Assemble dA
1675  try
1676  {
1677  tmpLocRhs = 0.0;
1678  Eval.add_def_A_elem(tmpLocRhs, locU, elem, vCornerCoords, PT_INSTATIONARY);
1679  locRhs.scale_append(-vScaleStiff[t], tmpLocRhs);
1680  }
1681  UG_CATCH_THROW("(instationary) AssembleLinear: Cannot compute Jacobian (A).");
1682 
1683  // Assemble rhs
1684  try
1685  {
1686  tmpLocRhs = 0.0;
1687  Eval.add_rhs_elem(tmpLocRhs, elem, vCornerCoords, PT_INSTATIONARY);
1688  locRhs.scale_append(vScaleStiff[t], tmpLocRhs);
1689  }
1690  UG_CATCH_THROW("(instationary) AssembleLinear: Cannot compute Rhs.");
1691  }
1692 
1693  // eliminate the outer connections in the local matrix
1694  m_extrapol.template eliminate_extrapolated<TElem> (locA, locRhs, base_co);
1695 
1696  // send local to global matrix & rhs
1697  try{
1698  spAssTuner->add_local_mat_to_global(A, locA, dd);
1699  spAssTuner->add_local_vec_to_global(rhs, locRhs, dd);
1700  }
1701  UG_CATCH_THROW("(instationary) AssembleLinear: Cannot add local vector/matrix.");
1702  }
1703  }
1704  }
1705 
1706 // finish element loop
1707  try
1708  {
1709  Eval.finish_elem_loop();
1710  }
1711  UG_CATCH_THROW("(instationary) AssembleLinear: Cannot finish element loop.");
1712 
1713  }
1714  UG_CATCH_THROW("(instationary) AssembleLinear: Cannot create Data Evaluator.");
1715 }
1716 
1717 template <typename TDomain, typename TAlgebra, typename TExtrapolation>
1720 (
1721  const std::vector<IElemDisc<domain_type>*>& vElemDisc,
1723  bool bNonRegularGrid,
1725  number future_time,
1727 )
1728 {
1729 // get current time and vector
1730  const vector_type& u = *vSol->solution(0);
1731 
1732 // create local time series
1733  LocalVectorTimeSeries locTimeSeries;
1734  locTimeSeries.read_times(vSol);
1735 
1736  try
1737  {
1739  vElemDisc, dd->function_pattern(), bNonRegularGrid,
1740  &locTimeSeries);
1741  Eval.set_time_point(0);
1742 
1743  // prepare timestep
1744  try
1745  {
1747  size_t algebra_index = bridge::AlgebraTypeIDProvider::instance().id<algebra_type>();
1748  Eval.prepare_timestep(future_time, vSol->time(0), &vp, algebra_index);
1749  }
1750  UG_CATCH_THROW("(instationary) PrepareTimestep: Cannot prepare timestep.");
1751 
1752  }
1753  UG_CATCH_THROW("(instationary) PrepareTimestep: Cannot create Data Evaluator.");
1754 }
1755 
1756 
1775 template <typename TDomain, typename TAlgebra, typename TExtrapolation>
1776 template <typename TElem, typename TIterator>
1778 PrepareTimestepElem(const std::vector<IElemDisc<domain_type>*>& vElemDisc,
1779  ConstSmartPtr<domain_type> spDomain,
1781  TIterator iterBegin,
1782  TIterator iterEnd,
1783  int si, bool bNonRegularGrid,
1786 {
1787 // check the dof distribution for the ghost-fluid method
1788  m_extrapol.check_dd (dd);
1789 
1790 // check if there are any elements at all, otherwise return immediately
1791  if(iterBegin == iterEnd) return;
1792 
1793 // reference object id
1795 
1796 // storage for corner coordinates
1798 
1799 // get current time and vector
1800  const vector_type& u = *vSol->solution(0);
1801 
1802 // create local time series
1803  LocalVectorTimeSeries locTimeSeries;
1804  locTimeSeries.read_times(vSol);
1805 
1806  try
1807  {
1809  vElemDisc, dd->function_pattern(), bNonRegularGrid,
1810  &locTimeSeries);
1811  Eval.set_time_point(0);
1812 
1813 // prepare element loop
1814  Eval.prepare_elem_loop(id, si);
1815 
1816 // local algebra
1817  LocalIndices ind; LocalVector locU;
1818 
1819 // grid level of all the elements
1820  int g_level = dd->grid_level().level ();
1821  if (g_level == GridLevel::TOP)
1822  g_level = dd->multi_grid()->top_level ();
1823 
1824 // Loop over all elements
1825  for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
1826  {
1827  // get Element
1828  TElem* elem = *iter;
1829 
1830  // check if elem is skipped from assembling
1831  if(!spAssTuner->element_used(elem)) continue;
1832 
1833  // get corner coordinates
1834  FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
1835 
1836  // check whether we are inside
1837  int elem_status = m_extrapol.check_elem_lsf
1838  (TElem::NUM_VERTICES, elem, si, g_level, Eval.use_hanging(), vCornerCoords, vSol->time(0));
1839  if (elem_status < 0) continue; // this is an outer element, do not assemble it at all
1840  if (m_bAssembleOnlyCut && elem_status > 0) continue; // exclude elements that are not cut
1841 
1842  // get global indices
1843  dd->indices(elem, ind, Eval.use_hanging());
1844 
1845  // adapt local algebra
1846  locU.resize(ind);
1847 
1848  // read local values of u
1849  GetLocalVector(locU, u);
1850 
1851  // read local values of time series
1852  if(Eval.time_series_needed())
1853  locTimeSeries.read_values(vSol, ind);
1854 
1855  // prepare timestep
1856  try
1857  {
1858  Eval.prepare_timestep_elem(vSol->time(0), locU, elem, vCornerCoords);
1859  }
1860  UG_CATCH_THROW("(instationary) PrepareTimestepElem: Cannot prepare timestep.");
1861  }
1862 
1863  }
1864  UG_CATCH_THROW("(instationary) PrepareTimestepElem: Cannot create Data Evaluator.");
1865 }
1866 
1867 template <typename TDomain, typename TAlgebra, typename TExtrapolation>
1870 (
1871  const std::vector<IElemDisc<domain_type>*>& vElemDisc,
1873  bool bNonRegularGrid,
1876 )
1877 {
1878 // get current time and vector
1879  const vector_type& u = *vSol->solution(0);
1880 
1881 // create local time series
1882  LocalVectorTimeSeries locTimeSeries;
1883  locTimeSeries.read_times(vSol);
1884 
1885  try
1886  {
1888  vElemDisc, dd->function_pattern(), bNonRegularGrid,
1889  &locTimeSeries);
1890  Eval.set_time_point(0);
1891 
1892  // finish time step
1893  try
1894  {
1896  size_t algebra_index = bridge::AlgebraTypeIDProvider::instance().id<algebra_type>();
1897  Eval.finish_timestep(vSol->time(0), &vp, algebra_index);
1898  }
1899  UG_CATCH_THROW("(instationary) FinishTimestep: Cannot prepare time step.");
1900 
1901  }
1902  UG_CATCH_THROW("(instationary) FinishTimestep: Cannot create Data Evaluator.");
1903 }
1904 
1922 template <typename TDomain, typename TAlgebra, typename TExtrapolation>
1923 template <typename TElem, typename TIterator>
1925 FinishTimestepElem(const std::vector<IElemDisc<domain_type>*>& vElemDisc,
1926  ConstSmartPtr<domain_type> spDomain,
1928  TIterator iterBegin,
1929  TIterator iterEnd,
1930  int si, bool bNonRegularGrid,
1933 {
1934 // check the dof distribution for the ghost-fluid method
1935  m_extrapol.check_dd (dd);
1936 
1937 // check if there are any elements at all, otherwise return immediately
1938  if(iterBegin == iterEnd) return;
1939 
1940 // reference object id
1942 
1943 // storage for corner coordinates
1945 
1946 // get current time and vector
1947  const vector_type& u = *vSol->solution(0);
1948 
1949 // create local time series
1950  LocalVectorTimeSeries locTimeSeries;
1951  locTimeSeries.read_times(vSol);
1952 
1953 // prepare for given elem discs
1954  try
1955  {
1957  vElemDisc, dd->function_pattern(), bNonRegularGrid,
1958  &locTimeSeries);
1959  Eval.set_time_point(0);
1960 
1961 // prepare loop
1962  Eval.prepare_elem_loop(id, si);
1963 
1964 // local algebra
1965  LocalIndices ind; LocalVector locU;
1966 
1967 // grid level of all the elements
1968  int g_level = dd->grid_level().level ();
1969  if (g_level == GridLevel::TOP)
1970  g_level = dd->multi_grid()->top_level ();
1971 
1972 // Loop over all elements
1973  for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
1974  {
1975  // get Element
1976  TElem* elem = *iter;
1977 
1978  // check if elem is skipped from assembling
1979  if(!spAssTuner->element_used(elem)) continue;
1980 
1981  // get corner coordinates
1982  FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
1983 
1984  // check whether we are inside
1985  int elem_status = m_extrapol.check_elem_lsf
1986  (TElem::NUM_VERTICES, elem, si, g_level, Eval.use_hanging(), vCornerCoords, vSol->time(0));
1987  if (elem_status < 0) continue; // this is an outer element, do not assemble it at all
1988  if (m_bAssembleOnlyCut && elem_status > 0) continue; // exclude elements that are not cut
1989 
1990  // get global indices
1991  dd->indices(elem, ind, Eval.use_hanging());
1992 
1993  // adapt local algebra
1994  locU.resize(ind);
1995 
1996  // read local values of u
1997  GetLocalVector(locU, u);
1998 
1999  // read local values of time series
2000  if(Eval.time_series_needed())
2001  locTimeSeries.read_values(vSol, ind);
2002 
2003  // finish timestep
2004  try
2005  {
2006  Eval.finish_timestep_elem(locTimeSeries.time(0), locU, elem, vCornerCoords);
2007  }
2008  UG_CATCH_THROW("(instationary) FinishTimestepElem: Cannot finish timestep.");
2009  }
2010 
2011  }
2012  UG_CATCH_THROW("(instationary) FinishTimestepElem: Cannot create Data Evaluator");
2013 }
2014 
2015 } // end namespace ug
2016 
2017 /* End of File */
The AssemblingTuner class combines tools to adapt the assembling routine.
Definition: ass_tuner.h:90
Evaluation for IElemDisc.
Definition: data_evaluator.h:169
@ TOP
Definition: grid_level.h:45
Definition: elem_disc_interface.h:800
void FinishTimestep(const std::vector< IElemDisc< domain_type > * > &vElemDisc, ConstSmartPtr< DoFDistribution > dd, bool bNonRegularGrid, ConstSmartPtr< VectorTimeSeries< vector_type > > vSol, ConstSmartPtr< AssemblingTuner< TAlgebra > > spAssTuner)
Definition: dom_disc_embb_impl.h:1870
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: dom_disc_embb_impl.h:1720
TAlgebra algebra_type
Algebra type.
Definition: dom_disc_embb.h:166
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: dom_disc_embb_impl.h:218
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: dom_disc_embb_impl.h:747
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: dom_disc_embb_impl.h:1264
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: dom_disc_embb_impl.h:60
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: dom_disc_embb_impl.h:1925
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: dom_disc_embb_impl.h:1778
algebra_type::vector_type vector_type
Vector type in the algebra.
Definition: dom_disc_embb.h:172
algebra_type::matrix_type matrix_type
Matrix type in the algebra.
Definition: dom_disc_embb.h:175
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: dom_disc_embb_impl.h:376
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
size_t size() const
returns number of time points
Definition: solution_time_series.h:173
a mathematical Vector with N entries.
Definition: math_vector.h:97
time series of solutions and corresponding time point
Definition: solution_time_series.h:59
Definition: grid_base_object_traits.h:68
#define UG_CATCH_THROW(msg)
Definition: error.h:64
#define UG_THROW(msg)
Definition: error.h:57
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
Definition: elem_disc_interface.h:77
size_t id()
Definition: util_algebra_dependent.h:130
static AlgebraTypeIDProvider & instance()
Definition: util_algebra_dependent.h:107