Line data Source code
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"
46 : #include "lib_disc/common/function_group.h"
47 : #include "lib_disc/common/local_algebra.h"
48 : #include "lib_disc/spatial_disc/user_data/data_evaluator.h"
49 : #include "bridge/util_algebra_dependent.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 :
65 : extern DebugID DID_ELEM_DISC_ASSEMBLE_UTIL;
66 :
67 : /// Global assembler based on the straightforward application of the local discretizations
68 : /**
69 : * Template class of the global assembler that applyes the local (element)
70 : * discretizations to the elements in given subsets and adds the local data
71 : * to the global ones.
72 : *
73 : * \tparam TDomain domain type
74 : * \tparam TAlgebra algebra type
75 : */
76 : template <typename TDomain, typename TAlgebra>
77 : class StdGlobAssembler
78 : {
79 : /// Domain type
80 : typedef TDomain domain_type;
81 :
82 : /// Algebra type
83 : typedef TAlgebra algebra_type;
84 :
85 : /// Vector type in the algebra
86 : typedef typename algebra_type::vector_type vector_type;
87 :
88 : /// Matrix type in the algebra
89 : typedef typename algebra_type::matrix_type matrix_type;
90 :
91 : ////////////////////////////////////////////////////////////////////////////////
92 : // Assemble Stiffness Matrix
93 : ////////////////////////////////////////////////////////////////////////////////
94 :
95 : public:
96 : /**
97 : * This function adds the contributions of all passed element discretizations
98 : * on one given subset to the global Stiffness matrix. (This version
99 : * processes elements in a given interval.)
100 : *
101 : * \param[in] vElemDisc element discretizations
102 : * \param[in] spDomain domain
103 : * \param[in] dd DoF Distribution
104 : * \param[in] iterBegin element iterator
105 : * \param[in] iterEnd element iterator
106 : * \param[in] si subset index
107 : * \param[in] bNonRegularGrid flag to indicate if non regular grid is used
108 : * \param[in,out] A Stiffness matrix
109 : * \param[in] u solution
110 : * \param[in] spAssTuner assemble adapter
111 : */
112 : template <typename TElem, typename TIterator>
113 : static void
114 0 : AssembleStiffnessMatrix( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
115 : ConstSmartPtr<domain_type> spDomain,
116 : ConstSmartPtr<DoFDistribution> dd,
117 : TIterator iterBegin,
118 : TIterator iterEnd,
119 : int si, bool bNonRegularGrid,
120 : matrix_type& A,
121 : const vector_type& u,
122 : ConstSmartPtr<AssemblingTuner<TAlgebra> > spAssTuner)
123 : {
124 : // check if there are any elements at all, otherwise return immediately
125 0 : if(iterBegin == iterEnd) return;
126 :
127 : // reference object id
128 : static const ReferenceObjectID id = geometry_traits<TElem>::REFERENCE_OBJECT_ID;
129 :
130 : // storage for corner coordinates
131 0 : MathVector<domain_type::dim> vCornerCoords[TElem::NUM_VERTICES];
132 :
133 : try
134 : {
135 0 : DataEvaluator<domain_type> Eval(STIFF,
136 : vElemDisc, dd->function_pattern(), bNonRegularGrid);
137 :
138 : // prepare element loop
139 0 : 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 0 : for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
146 : {
147 : // get Element
148 0 : TElem* elem = *iter;
149 :
150 : // get corner coordinates
151 0 : FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
152 :
153 : // check if elem is skipped from assembling
154 0 : if(!spAssTuner->element_used(elem)) continue;
155 :
156 : // get global indices
157 0 : dd->indices(elem, ind, Eval.use_hanging());
158 :
159 : // adapt local algebra
160 0 : locU.resize(ind); locA.resize(ind);
161 :
162 : // read local values of u
163 0 : GetLocalVector(locU, u);
164 :
165 : // prepare element
166 : try
167 : {
168 0 : Eval.prepare_elem(locU, elem, id, vCornerCoords, ind, true);
169 : }
170 0 : UG_CATCH_THROW("AssembleStiffnessMatrix Cannot prepare element.");
171 :
172 : // Assemble JA
173 0 : locA = 0.0;
174 : try
175 : {
176 0 : Eval.add_jac_A_elem(locA, locU, elem, vCornerCoords);
177 : }
178 0 : UG_CATCH_THROW("AssembleStiffnessMatrix: Cannot compute Jacobian (A).");
179 :
180 : // send local to global matrix
181 : try{
182 0 : spAssTuner->add_local_mat_to_global(A,locA,dd);
183 : }
184 0 : UG_CATCH_THROW("AssembleStiffnessMatrix: Cannot add local matrix.");
185 : }
186 :
187 : // finish element loop
188 : try
189 : {
190 0 : Eval.finish_elem_loop();
191 : }
192 0 : UG_CATCH_THROW("AssembleStiffnessMatrix: Cannot finish element loop.");
193 :
194 0 : }
195 0 : UG_CATCH_THROW("AssembleStiffnessMatrix': Cannot create Data Evaluator.");
196 : }
197 :
198 : ////////////////////////////////////////////////////////////////////////////////
199 : // Assemble Mass Matrix
200 : ////////////////////////////////////////////////////////////////////////////////
201 :
202 : public:
203 : /**
204 : * This function adds the contributions of all passed element discretizations
205 : * on one given subset to the global Mass matrix. (This version
206 : * processes elements in a given interval.)
207 : *
208 : * \param[in] vElemDisc element discretizations
209 : * \param[in] spDomain domain
210 : * \param[in] dd DoF Distribution
211 : * \param[in] iterBegin element iterator
212 : * \param[in] iterEnd element iterator
213 : * \param[in] si subset index
214 : * \param[in] bNonRegularGrid flag to indicate if non regular grid is used
215 : * \param[in,out] M Mass matrix
216 : * \param[in] u solution
217 : * \param[in] spAssTuner assemble adapter
218 : */
219 : template <typename TElem, typename TIterator>
220 : static void
221 0 : AssembleMassMatrix( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
222 : ConstSmartPtr<domain_type> spDomain,
223 : ConstSmartPtr<DoFDistribution> dd,
224 : TIterator iterBegin,
225 : TIterator iterEnd,
226 : int si, bool bNonRegularGrid,
227 : matrix_type& M,
228 : const vector_type& u,
229 : ConstSmartPtr<AssemblingTuner<TAlgebra> > spAssTuner)
230 : {
231 : // check if there are any elements at all, otherwise return immediately
232 0 : if(iterBegin == iterEnd) return;
233 :
234 : // reference object id
235 : static const ReferenceObjectID id = geometry_traits<TElem>::REFERENCE_OBJECT_ID;
236 :
237 : // storage for corner coordinates
238 0 : MathVector<domain_type::dim> vCornerCoords[TElem::NUM_VERTICES];
239 :
240 : // prepare for given elem discs
241 : try
242 : {
243 0 : DataEvaluator<domain_type> Eval(MASS,
244 : vElemDisc, dd->function_pattern(), bNonRegularGrid);
245 :
246 : // prepare element loop
247 0 : 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 0 : for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
254 : {
255 : // get Element
256 0 : TElem* elem = *iter;
257 :
258 : // get corner coordinates
259 0 : FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
260 :
261 : // check if elem is skipped from assembling
262 0 : if(!spAssTuner->element_used(elem)) continue;
263 :
264 : // get global indices
265 0 : dd->indices(elem, ind, Eval.use_hanging());
266 :
267 : // adapt local algebra
268 0 : locU.resize(ind); locM.resize(ind);
269 :
270 : // read local values of u
271 0 : GetLocalVector(locU, u);
272 :
273 : // prepare element
274 : try
275 : {
276 0 : Eval.prepare_elem(locU, elem, id, vCornerCoords, ind, true);
277 : }
278 0 : UG_CATCH_THROW("AssembleMassMatrix: Cannot prepare element.");
279 :
280 : // Assemble JM
281 0 : locM = 0.0;
282 : try
283 : {
284 0 : Eval.add_jac_M_elem(locM, locU, elem, vCornerCoords);
285 : }
286 0 : UG_CATCH_THROW("AssembleMassMatrix: Cannot compute Jacobian (M).");
287 :
288 : // send local to global matrix
289 : try{
290 0 : spAssTuner->add_local_mat_to_global(M, locM, dd);
291 : }
292 0 : UG_CATCH_THROW("AssembleMassMatrix: Cannot add local matrix.");
293 : }
294 :
295 : // finish element loop
296 : try
297 : {
298 0 : Eval.finish_elem_loop();
299 : }
300 0 : UG_CATCH_THROW("AssembleMassMatrix: Cannot finish element loop.");
301 :
302 0 : }
303 0 : UG_CATCH_THROW("AssembleMassMatrix: Cannot create Data Evaluator.");
304 : }
305 :
306 : ////////////////////////////////////////////////////////////////////////////////
307 : // Assemble (stationary) Jacobian
308 : ////////////////////////////////////////////////////////////////////////////////
309 :
310 : public:
311 : /**
312 : * This function adds the contributions of all passed element discretizations
313 : * on one given subset to the global Jacobian in the stationary case. (This version
314 : * processes elements in a given interval.)
315 : *
316 : * \param[in] vElemDisc element discretizations
317 : * \param[in] spDomain domain
318 : * \param[in] iterBegin element iterator
319 : * \param[in] iterEnd element iterator
320 : * \param[in] si subset index
321 : * \param[in] bNonRegularGrid flag to indicate if non regular grid is used
322 : * \param[in,out] J jacobian
323 : * \param[in] u solution
324 : * \param[in] spAssTuner assemble adapter
325 : */
326 : template <typename TElem, typename TIterator>
327 : static void
328 0 : AssembleJacobian( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
329 : ConstSmartPtr<domain_type> spDomain,
330 : ConstSmartPtr<DoFDistribution> dd,
331 : TIterator iterBegin,
332 : TIterator iterEnd,
333 : int si, bool bNonRegularGrid,
334 : matrix_type& J,
335 : const vector_type& u,
336 : ConstSmartPtr<AssemblingTuner<TAlgebra> > spAssTuner)
337 : {
338 : // check if there are any elements at all, otherwise return immediately
339 0 : if(iterBegin == iterEnd) return;
340 :
341 : // reference object id
342 : static const ReferenceObjectID id = geometry_traits<TElem>::REFERENCE_OBJECT_ID;
343 :
344 : // storage for corner coordinates
345 0 : MathVector<domain_type::dim> vCornerCoords[TElem::NUM_VERTICES];
346 :
347 : // prepare for given elem discs
348 : try
349 : {
350 0 : DataEvaluator<domain_type> Eval(STIFF | RHS,
351 : vElemDisc, dd->function_pattern(), bNonRegularGrid);
352 :
353 : // prepare element loop
354 0 : 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 0 : for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
361 : {
362 : // get Element
363 0 : TElem* elem = *iter;
364 :
365 : // get corner coordinates
366 0 : FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
367 :
368 : // check if elem is skipped from assembling
369 0 : if(!spAssTuner->element_used(elem)) continue;
370 :
371 : // get global indices
372 0 : dd->indices(elem, ind, Eval.use_hanging());
373 :
374 : // adapt local algebra
375 0 : locU.resize(ind); locJ.resize(ind);
376 :
377 : // read local values of u
378 0 : GetLocalVector(locU, u);
379 :
380 : // prepare element
381 : try
382 : {
383 0 : Eval.prepare_elem(locU, elem, id, vCornerCoords, ind, true);
384 : }
385 0 : UG_CATCH_THROW("(stationary) AssembleJacobian: Cannot prepare element.");
386 :
387 : // reset local algebra
388 0 : locJ = 0.0;
389 :
390 : // Assemble JA
391 : try
392 : {
393 0 : Eval.add_jac_A_elem(locJ, locU, elem, vCornerCoords);
394 : }
395 0 : UG_CATCH_THROW("(stationary) AssembleJacobian: Cannot compute Jacobian (A).");
396 :
397 : // send local to global matrix
398 : try{
399 0 : spAssTuner->add_local_mat_to_global(J, locJ, dd);
400 : }
401 0 : UG_CATCH_THROW("(stationary) AssembleJacobian: Cannot add local matrix.");
402 : }
403 :
404 : // finish element loop
405 : try
406 : {
407 0 : Eval.finish_elem_loop();
408 : }
409 0 : UG_CATCH_THROW("(stationary) AssembleJacobian: Cannot finish element loop.");
410 :
411 0 : }
412 0 : UG_CATCH_THROW("(stationary) AssembleJacobian: Cannot create Data Evaluator.");
413 : }
414 :
415 : ////////////////////////////////////////////////////////////////////////////////
416 : // Assemble (instationary) Jacobian
417 : ////////////////////////////////////////////////////////////////////////////////
418 :
419 : public:
420 : /**
421 : * This function adds the contributions of all passed element discretizations
422 : * on one given subset to the global Jacobian in the time-dependent case.
423 : * Note, that s_m0 == 1
424 : * (This version processes elements in a given interval.)
425 : *
426 : * \param[in] vElemDisc element discretizations
427 : * \param[in] spDomain domain
428 : * \param[in] dd DoF Distribution
429 : * \param[in] iterBegin element iterator
430 : * \param[in] iterEnd element iterator
431 : * \param[in] si subset index
432 : * \param[in] bNonRegularGrid flag to indicate if non regular grid is used
433 : * \param[in,out] J jacobian
434 : * \param[in] vSol current and previous solutions
435 : * \param[in] s_a0 scaling factor for stiffness part
436 : * \param[in] spAssTuner assemble adapter
437 : */
438 : template <typename TElem, typename TIterator>
439 : static void
440 0 : AssembleJacobian( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
441 : ConstSmartPtr<domain_type> spDomain,
442 : ConstSmartPtr<DoFDistribution> dd,
443 : TIterator iterBegin,
444 : TIterator iterEnd,
445 : int si, bool bNonRegularGrid,
446 : matrix_type& J,
447 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
448 : number s_a0,
449 : ConstSmartPtr<AssemblingTuner<TAlgebra> > spAssTuner)
450 : {
451 : // check if there are any elements at all, otherwise return immediately
452 0 : if(iterBegin == iterEnd) return;
453 :
454 : // reference object id
455 : static const ReferenceObjectID id = geometry_traits<TElem>::REFERENCE_OBJECT_ID;
456 :
457 : // storage for corner coordinates
458 0 : MathVector<domain_type::dim> vCornerCoords[TElem::NUM_VERTICES];
459 :
460 : // get current time and vector
461 0 : const vector_type& u = *vSol->solution(0);
462 :
463 : // create local time series
464 : LocalVectorTimeSeries locTimeSeries;
465 0 : locTimeSeries.read_times(vSol);
466 :
467 : // prepare for given elem discs
468 : try
469 : {
470 0 : DataEvaluator<domain_type> Eval(MASS | STIFF | RHS,
471 : vElemDisc, dd->function_pattern(), bNonRegularGrid,
472 : &locTimeSeries);
473 0 : Eval.set_time_point(0);
474 :
475 : // prepare element loop
476 0 : 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 0 : for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
484 : {
485 : // get Element
486 0 : TElem* elem = *iter;
487 :
488 : // get corner coordinates
489 0 : FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
490 :
491 : // check if elem is skipped from assembling
492 0 : if(!spAssTuner->element_used(elem)) continue;
493 :
494 : // get global indices
495 0 : dd->indices(elem, ind, Eval.use_hanging());
496 :
497 : // adapt local algebra
498 0 : locU.resize(ind); locJ.resize(ind);
499 :
500 : // read local values of u
501 0 : GetLocalVector(locU, u);
502 :
503 : // read local values of time series
504 0 : if(Eval.time_series_needed())
505 0 : locTimeSeries.read_values(vSol, ind);
506 :
507 : // prepare element
508 : try
509 : {
510 0 : Eval.prepare_elem(locU, elem, id, vCornerCoords, ind, true);
511 : }
512 0 : UG_CATCH_THROW("(instationary) AssembleJacobian: Cannot prepare element.");
513 :
514 : // reset local algebra
515 0 : locJ = 0.0;
516 :
517 : //EL_PROFILE_BEGIN(Elem_add_JA);
518 : // Assemble JA
519 : try
520 : {
521 0 : Eval.add_jac_A_elem(locJ, locU, elem, vCornerCoords, PT_INSTATIONARY);
522 0 : locJ *= s_a0;
523 :
524 0 : Eval.add_jac_A_elem(locJ, locU, elem, vCornerCoords, PT_STATIONARY);
525 : }
526 0 : UG_CATCH_THROW("(instationary) AssembleJacobian: Cannot compute Jacobian (A).");
527 : //EL_PROFILE_END();
528 :
529 : // Assemble JM
530 : try
531 : {
532 0 : Eval.add_jac_M_elem(locJ, locU, elem, vCornerCoords, PT_INSTATIONARY);
533 : }
534 0 : UG_CATCH_THROW("(instationary) AssembleJacobian: Cannot compute Jacobian (M).");
535 :
536 : // send local to global matrix
537 : try{
538 0 : spAssTuner->add_local_mat_to_global(J, locJ, dd);
539 : }
540 0 : UG_CATCH_THROW("(instationary) AssembleJacobian: Cannot add local matrix.");
541 :
542 : }
543 : EL_PROFILE_END();
544 :
545 : // finish element loop
546 : try
547 : {
548 0 : Eval.finish_elem_loop();
549 : }
550 0 : UG_CATCH_THROW("(instationary) AssembleJacobian: Cannot finish element loop.");
551 :
552 0 : }
553 0 : UG_CATCH_THROW("(instationary) AssembleJacobian: Cannot create Data Evaluator.");
554 : }
555 :
556 : ////////////////////////////////////////////////////////////////////////////////
557 : // Assemble (stationary) Defect
558 : ////////////////////////////////////////////////////////////////////////////////
559 :
560 : public:
561 : /**
562 : * This function adds the contributions of all passed element discretizations
563 : * on one given subset to the global Defect in the stationary case. (This version
564 : * processes elements in a given interval.)
565 : *
566 : * \param[in] vElemDisc element discretizations
567 : * \param[in] spDomain domain
568 : * \param[in] dd DoF Distribution
569 : * \param[in] iterBegin element iterator
570 : * \param[in] iterEnd element iterator
571 : * \param[in] si subset index
572 : * \param[in] bNonRegularGrid flag to indicate if non regular grid is used
573 : * \param[in,out] d defect
574 : * \param[in] u solution
575 : * \param[in] spAssTuner assemble adapter
576 : */
577 : template <typename TElem, typename TIterator>
578 : static void
579 0 : AssembleDefect( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
580 : ConstSmartPtr<domain_type> spDomain,
581 : ConstSmartPtr<DoFDistribution> dd,
582 : TIterator iterBegin,
583 : TIterator iterEnd,
584 : int si, bool bNonRegularGrid,
585 : vector_type& d,
586 : const vector_type& u,
587 : ConstSmartPtr<AssemblingTuner<TAlgebra> > spAssTuner)
588 : {
589 : // check if at least one element exists, else return
590 0 : if(iterBegin == iterEnd) return;
591 :
592 : // reference object id
593 : static const ReferenceObjectID id = geometry_traits<TElem>::REFERENCE_OBJECT_ID;
594 :
595 : // storage for corner coordinates
596 0 : MathVector<domain_type::dim> vCornerCoords[TElem::NUM_VERTICES];
597 :
598 : // prepare for given elem discs
599 : try
600 : {
601 0 : DataEvaluator<domain_type> Eval(STIFF | RHS,
602 : vElemDisc, dd->function_pattern(), bNonRegularGrid);
603 :
604 : // prepare element loop
605 0 : 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 0 : for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
612 : {
613 : // get Element
614 0 : TElem* elem = *iter;
615 :
616 : // get corner coordinates
617 0 : FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
618 :
619 : // check if elem is skipped from assembling
620 0 : if(!spAssTuner->element_used(elem)) continue;
621 :
622 : // get global indices
623 0 : dd->indices(elem, ind, Eval.use_hanging());
624 :
625 : // adapt local algebra
626 0 : locU.resize(ind); locD.resize(ind); tmpLocD.resize(ind);
627 :
628 : // read local values of u
629 0 : GetLocalVector(locU, u);
630 :
631 : // prepare element
632 : try
633 : {
634 0 : Eval.prepare_elem(locU, elem, id, vCornerCoords, ind);
635 : }
636 0 : 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 0 : if( spAssTuner->modify_solution_enabled() )
641 : {
642 : LocalVector& modLocU = locU;
643 : try{
644 0 : spAssTuner->modify_LocalSol(modLocU, locU, dd);
645 0 : } UG_CATCH_THROW("Cannot modify local solution.");
646 :
647 : // recopy modified LocalVector:
648 0 : locU = modLocU;
649 : }
650 :
651 : // reset local algebra
652 0 : locD = 0.0;
653 :
654 : // Assemble A
655 : try
656 : {
657 0 : Eval.add_def_A_elem(locD, locU, elem, vCornerCoords);
658 : }
659 0 : UG_CATCH_THROW("(stationary) AssembleDefect: Cannot compute Defect (A).");
660 :
661 : // Assemble rhs
662 : try
663 : {
664 0 : tmpLocD = 0.0;
665 0 : Eval.add_rhs_elem(tmpLocD, elem, vCornerCoords);
666 0 : locD.scale_append(-1, tmpLocD);
667 :
668 : }
669 0 : UG_CATCH_THROW("(stationary) AssembleDefect: Cannot compute Rhs.");
670 :
671 : // send local to global defect
672 : try{
673 0 : spAssTuner->add_local_vec_to_global(d, locD, dd);
674 : }
675 0 : UG_CATCH_THROW("(stationary) AssembleDefect: Cannot add local vector.");
676 : }
677 :
678 : // finish element loop
679 : try
680 : {
681 0 : Eval.finish_elem_loop();
682 : }
683 0 : UG_CATCH_THROW("(stationary) AssembleDefect: Cannot finish element loop.");
684 :
685 : }
686 0 : UG_CATCH_THROW("(stationary) AssembleDefect: Cannot create Data Evaluator.");
687 : }
688 :
689 : ////////////////////////////////////////////////////////////////////////////////
690 : // Assemble (instationary) Defect
691 : ////////////////////////////////////////////////////////////////////////////////
692 :
693 : public:
694 : /**
695 : * This function adds the contributions of all passed element discretizations
696 : * on one given subset to the global Defect in the instationary case. (This version
697 : * processes elements in a given interval.)
698 : *
699 : * \param[in] vElemDisc element discretizations
700 : * \param[in] spDomain domain
701 : * \param[in] dd DoF Distribution
702 : * \param[in] iterBegin element iterator
703 : * \param[in] iterEnd element iterator
704 : * \param[in] si subset index
705 : * \param[in] bNonRegularGrid flag to indicate if non regular grid is used
706 : * \param[in,out] d defect
707 : * \param[in] vSol current and previous solutions
708 : * \param[in] vScaleMass scaling factors for mass part
709 : * \param[in] vScaleStiff scaling factors for stiffness part
710 : * \param[in] spAssTuner assemble adapter
711 : */
712 : template <typename TElem, typename TIterator>
713 : static void
714 0 : AssembleDefect( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
715 : ConstSmartPtr<domain_type> spDomain,
716 : ConstSmartPtr<DoFDistribution> dd,
717 : TIterator iterBegin,
718 : TIterator iterEnd,
719 : int si, bool bNonRegularGrid,
720 : vector_type& d,
721 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
722 : const std::vector<number>& vScaleMass,
723 : const std::vector<number>& vScaleStiff,
724 : ConstSmartPtr<AssemblingTuner<TAlgebra> > spAssTuner)
725 : {
726 : // check if there are any elements at all, otherwise return immediately
727 0 : if(iterBegin == iterEnd) return;
728 :
729 : // reference object id
730 : static const ReferenceObjectID id = geometry_traits<TElem>::REFERENCE_OBJECT_ID;
731 :
732 : // storage for corner coordinates
733 0 : MathVector<domain_type::dim> vCornerCoords[TElem::NUM_VERTICES];
734 :
735 : // check time scheme
736 0 : if(vScaleMass.size() != vScaleStiff.size())
737 0 : UG_THROW("(instationary) AssembleDefect: s_a and s_m must have same size.");
738 :
739 0 : if(vSol->size() < vScaleStiff.size())
740 0 : 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 0 : locTimeSeries.read_times(vSol);
747 :
748 : // prepare for given elem discs
749 : try
750 : {
751 0 : DataEvaluator<domain_type> Eval(MASS | STIFF | RHS | EXPL,
752 : vElemDisc, dd->function_pattern(), bNonRegularGrid,
753 : &locTimeSeries, &vScaleMass, &vScaleStiff);
754 :
755 : // prepare element loop
756 0 : 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 0 : for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
763 : {
764 : // get Element
765 0 : TElem* elem = *iter;
766 :
767 : // get corner coordinates
768 0 : FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
769 :
770 : // check if elem is skipped from assembling
771 0 : if(!spAssTuner->element_used(elem)) continue;
772 :
773 : // get global indices
774 0 : dd->indices(elem, ind, Eval.use_hanging());
775 :
776 : // adapt local algebra
777 0 : locD.resize(ind); tmpLocD.resize(ind);
778 :
779 : // read local values of time series
780 0 : locTimeSeries.read_values(vSol, ind);
781 :
782 : // reset contribution of this element
783 0 : locD = 0.0;
784 :
785 : // loop all time points and assemble them
786 0 : for(size_t t = 0; t < vScaleStiff.size(); ++t)
787 : {
788 0 : number scale_stiff = vScaleStiff[t];
789 :
790 : // get local solution at timepoint
791 : LocalVector& locU = locTimeSeries.solution(t);
792 0 : Eval.set_time_point(t);
793 :
794 : // prepare element
795 : try
796 : {
797 0 : Eval.prepare_elem(locU, elem, id, vCornerCoords, ind, false);
798 : }
799 0 : UG_CATCH_THROW("(instationary) AssembleDefect: Cannot prepare element.");
800 :
801 : // Assemble M
802 : try
803 : {
804 0 : tmpLocD = 0.0;
805 0 : Eval.add_def_M_elem(tmpLocD, locU, elem, vCornerCoords, PT_INSTATIONARY);
806 0 : locD.scale_append(vScaleMass[t], tmpLocD);
807 : }
808 0 : UG_CATCH_THROW("(instationary) AssembleDefect: Cannot compute Defect (M).");
809 :
810 : // Assemble A
811 : try
812 : {
813 0 : if(scale_stiff != 0.0)
814 : {
815 0 : tmpLocD = 0.0;
816 0 : Eval.add_def_A_elem(tmpLocD, locU, elem, vCornerCoords, PT_INSTATIONARY);
817 0 : locD.scale_append(scale_stiff, tmpLocD);
818 : }
819 :
820 0 : if(t == 0)
821 0 : Eval.add_def_A_elem(locD, locU, elem, vCornerCoords, PT_STATIONARY);
822 : }
823 0 : UG_CATCH_THROW("(instationary) AssembleDefect: Cannot compute Defect (A).");
824 :
825 :
826 : // Assemble defect for explicit reaction_rate, reaction and source
827 0 : if( t == 1 ) // only valid at lowest timediscretization order
828 : {
829 0 : tmpLocD = 0.0;
830 : try
831 : {
832 0 : Eval.add_def_A_expl_elem(tmpLocD, locU, elem, vCornerCoords, PT_INSTATIONARY);
833 : }
834 0 : UG_CATCH_THROW("(instationary) AssembleDefect explizit: Cannot compute Defect (A).");
835 :
836 : // UG_ASSERT(vScaleStiff.size() == 2, "Only one step method supported.");
837 0 : const number dt = vSol->time(0)-vSol->time(1);
838 0 : locD.scale_append(dt, tmpLocD);
839 : }
840 :
841 :
842 : // Assemble rhs
843 : try
844 : {
845 0 : if(scale_stiff != 0.0)
846 : {
847 0 : tmpLocD = 0.0;
848 0 : Eval.add_rhs_elem(tmpLocD, elem, vCornerCoords, PT_INSTATIONARY);
849 0 : locD.scale_append( - scale_stiff, tmpLocD);
850 : }
851 :
852 0 : if(t == 0)
853 : {
854 0 : tmpLocD = 0.0;
855 0 : Eval.add_rhs_elem(tmpLocD, elem, vCornerCoords, PT_STATIONARY);
856 0 : locD.scale_append( -1.0, tmpLocD);
857 : }
858 : }
859 0 : UG_CATCH_THROW("(instationary) AssembleDefect: Cannot compute Rhs.");
860 : }
861 :
862 : // send local to global defect
863 : try{
864 0 : spAssTuner->add_local_vec_to_global(d, locD, dd);
865 : }
866 0 : UG_CATCH_THROW("(instationary) AssembleDefect: Cannot add local vector.");
867 : }
868 :
869 : // finish element loop
870 : try
871 : {
872 0 : Eval.finish_elem_loop();
873 : }
874 0 : UG_CATCH_THROW("(instationary) AssembleDefect: Cannot finish element loop.");
875 :
876 : }
877 0 : UG_CATCH_THROW("(instationary) AssembleDefect: Cannot create Data Evaluator.");
878 : }
879 :
880 : ////////////////////////////////////////////////////////////////////////////////
881 : // Assemble (stationary) Linear problem
882 : ////////////////////////////////////////////////////////////////////////////////
883 :
884 : public:
885 : /**
886 : * This function adds the contributions of all passed element discretizations
887 : * on one given subset to the global Matrix and the global Right-Hand Side
888 : * of the Linear problem in the stationary case. (This version processes
889 : * elements in a given interval.)
890 : *
891 : * \param[in] vElemDisc element discretizations
892 : * \param[in] spDomain domain
893 : * \param[in] dd DoF Distribution
894 : * \param[in] iterBegin element iterator
895 : * \param[in] iterEnd element iterator
896 : * \param[in] si subset index
897 : * \param[in] bNonRegularGrid flag to indicate if non regular grid is used
898 : * \param[in,out] A Matrix
899 : * \param[in,out] rhs Right-hand side
900 : * \param[in] spAssTuner assemble adapter
901 : */
902 : template <typename TElem, typename TIterator>
903 : static void
904 0 : AssembleLinear( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
905 : ConstSmartPtr<domain_type> spDomain,
906 : ConstSmartPtr<DoFDistribution> dd,
907 : TIterator iterBegin,
908 : TIterator iterEnd,
909 : int si, bool bNonRegularGrid,
910 : matrix_type& A,
911 : vector_type& rhs,
912 : ConstSmartPtr<AssemblingTuner<TAlgebra> > spAssTuner)
913 : {
914 : // check if there are any elements at all, otherwise return immediately
915 0 : if(iterBegin == iterEnd) return;
916 :
917 : // reference object id
918 : static const ReferenceObjectID id = geometry_traits<TElem>::REFERENCE_OBJECT_ID;
919 :
920 : // storage for corner coordinates
921 0 : MathVector<domain_type::dim> vCornerCoords[TElem::NUM_VERTICES];
922 :
923 : // prepare for given elem discs
924 : try
925 : {
926 0 : DataEvaluator<domain_type> Eval(STIFF | RHS,
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 0 : 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 0 : for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
941 : {
942 : // get Element
943 0 : TElem* elem = *iter;
944 :
945 : // get corner coordinates
946 0 : FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
947 :
948 : // check if elem is skipped from assembling
949 0 : if(!spAssTuner->element_used(elem)) continue;
950 :
951 : // get global indices
952 0 : dd->indices(elem, ind, Eval.use_hanging());
953 :
954 : // adapt local algebra
955 0 : 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 0 : Eval.prepare_elem(locRhs, elem, id, vCornerCoords, ind, true);
965 : }
966 0 : UG_CATCH_THROW("(stationary) AssembleLinear: Cannot prepare element.");
967 :
968 : // reset local algebra
969 0 : locA = 0.0;
970 0 : locRhs = 0.0;
971 :
972 : // Assemble JA
973 : try
974 : {
975 0 : Eval.add_jac_A_elem(locA, locRhs, elem, vCornerCoords);
976 : }
977 0 : UG_CATCH_THROW("(stationary) AssembleLinear: Cannot compute Jacobian (A).");
978 :
979 : // Assemble rhs
980 : try
981 : {
982 0 : Eval.add_rhs_elem(locRhs, elem, vCornerCoords);
983 : }
984 0 : UG_CATCH_THROW("(stationary) AssembleLinear: Cannot compute Rhs.");
985 :
986 : // send local to global matrix & rhs
987 : try{
988 0 : spAssTuner->add_local_mat_to_global(A, locA, dd);
989 0 : spAssTuner->add_local_vec_to_global(rhs, locRhs, dd);
990 : }
991 0 : UG_CATCH_THROW("(stationary) AssembleLinear: Cannot add local vector/matrix.");
992 : }
993 :
994 : // finish element loop
995 : try
996 : {
997 0 : Eval.finish_elem_loop();
998 : }
999 0 : UG_CATCH_THROW("(stationary) AssembleLinear: Cannot finish element loop.");
1000 :
1001 0 : }
1002 0 : UG_CATCH_THROW("(stationary) AssembleLinear: Cannot create Data Evaluator.");
1003 : }
1004 :
1005 : ////////////////////////////////////////////////////////////////////////////////
1006 : // Assemble (instationary) Linear problem
1007 : ////////////////////////////////////////////////////////////////////////////////
1008 :
1009 : public:
1010 : /**
1011 : * This function adds the contributions of all passed element discretizations
1012 : * on one given subset to the global Matrix and the global Right-Hand Side
1013 : * of the Linear problem in the time-dependent case. (This version processes
1014 : * elements in a given interval.)
1015 : *
1016 : * \param[in] vElemDisc element discretizations
1017 : * \param[in] spDomain domain
1018 : * \param[in] dd DoF Distribution
1019 : * \param[in] iterBegin element iterator
1020 : * \param[in] iterEnd element iterator
1021 : * \param[in] si subset index
1022 : * \param[in] bNonRegularGrid flag to indicate if non regular grid is used
1023 : * \param[in,out] A Matrix
1024 : * \param[in,out] rhs Right-hand side
1025 : * \param[in] vSol current and previous solutions
1026 : * \param[in] vScaleMass scaling factors for mass part
1027 : * \param[in] vScaleStiff scaling factors for stiffness part
1028 : * \param[in] spAssTuner assemble adapter
1029 : */
1030 : template <typename TElem, typename TIterator>
1031 : static void
1032 0 : AssembleLinear( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
1033 : ConstSmartPtr<domain_type> spDomain,
1034 : ConstSmartPtr<DoFDistribution> dd,
1035 : TIterator iterBegin,
1036 : TIterator iterEnd,
1037 : int si, bool bNonRegularGrid,
1038 : matrix_type& A,
1039 : vector_type& rhs,
1040 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
1041 : const std::vector<number>& vScaleMass,
1042 : const std::vector<number>& vScaleStiff,
1043 : ConstSmartPtr<AssemblingTuner<TAlgebra> > spAssTuner)
1044 : {
1045 : // check if there are any elements at all, otherwise return immediately
1046 0 : if(iterBegin == iterEnd) return;
1047 :
1048 : // reference object id
1049 : static const ReferenceObjectID id = geometry_traits<TElem>::REFERENCE_OBJECT_ID;
1050 :
1051 : // storage for corner coordinates
1052 0 : MathVector<domain_type::dim> vCornerCoords[TElem::NUM_VERTICES];
1053 :
1054 : // check time scheme
1055 0 : if(vScaleMass.size() != vScaleStiff.size())
1056 0 : UG_THROW("(instationary) AssembleLinear: s_a and s_m must have same size.");
1057 :
1058 0 : if(vSol->size() < vScaleStiff.size())
1059 0 : 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 0 : locTimeSeries.read_times(vSol);
1066 :
1067 : // prepare for given elem discs
1068 : try
1069 : {
1070 0 : DataEvaluator<domain_type> Eval(MASS | STIFF | RHS,
1071 : vElemDisc, dd->function_pattern(), bNonRegularGrid,
1072 : &locTimeSeries, &vScaleMass, &vScaleStiff);
1073 :
1074 : // prepare loop
1075 0 : 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 0 : for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
1082 : {
1083 : // get Element
1084 0 : TElem* elem = *iter;
1085 :
1086 : // get corner coordinates
1087 0 : FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
1088 :
1089 : // check if elem is skipped from assembling
1090 0 : if(!spAssTuner->element_used(elem)) continue;
1091 :
1092 : // get global indices
1093 0 : dd->indices(elem, ind, Eval.use_hanging());
1094 :
1095 : // adapt local algebra
1096 0 : locRhs.resize(ind); tmpLocRhs.resize(ind);
1097 : locA.resize(ind); tmpLocA.resize(ind);
1098 :
1099 : // read local values of time series
1100 0 : locTimeSeries.read_values(vSol, ind);
1101 0 : Eval.set_time_point(0);
1102 :
1103 : // reset element contribution
1104 0 : locA = 0.0; locRhs = 0.0;
1105 :
1106 : /////////////////////
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 0 : Eval.prepare_elem(locU, elem, id, vCornerCoords, ind, true);
1116 : }
1117 0 : UG_CATCH_THROW("(instationary) AssembleLinear: Cannot prepare element.");
1118 :
1119 0 : if (!spAssTuner->matrix_is_const())
1120 : {
1121 : // Assemble JM
1122 : try
1123 : {
1124 0 : tmpLocA = 0.0;
1125 0 : Eval.add_jac_M_elem(tmpLocA, locU, elem, vCornerCoords, PT_INSTATIONARY);
1126 0 : locA.scale_append(vScaleMass[0], tmpLocA);
1127 : }
1128 0 : UG_CATCH_THROW("(instationary) AssembleLinear: Cannot compute Jacobian (M).");
1129 :
1130 : // Assemble JA
1131 : try
1132 : {
1133 0 : if (vScaleStiff[0] != 0.0)
1134 : {
1135 0 : tmpLocA = 0.0;
1136 0 : Eval.add_jac_A_elem(tmpLocA, locU, elem, vCornerCoords, PT_INSTATIONARY);
1137 0 : locA.scale_append(vScaleStiff[0], tmpLocA);
1138 : }
1139 0 : Eval.add_jac_A_elem(locA, locU, elem, vCornerCoords, PT_STATIONARY);
1140 : }
1141 0 : UG_CATCH_THROW("(instationary) AssembleLinear: Cannot compute Jacobian (A).");
1142 : }
1143 : // Assemble rhs
1144 : try
1145 : {
1146 0 : if (vScaleStiff[0] != 0.0)
1147 : {
1148 0 : tmpLocRhs = 0.0;
1149 0 : Eval.add_rhs_elem(tmpLocRhs, elem, vCornerCoords, PT_INSTATIONARY);
1150 0 : locRhs.scale_append(vScaleStiff[0], tmpLocRhs);
1151 : }
1152 0 : Eval.add_rhs_elem(locRhs, elem, vCornerCoords, PT_STATIONARY);
1153 : }
1154 0 : UG_CATCH_THROW("(instationary) AssembleLinear: Cannot compute Rhs.");
1155 :
1156 : ///////////////////
1157 : // old time steps
1158 :
1159 : // loop all old time points
1160 0 : for(size_t t = 1; t < vScaleStiff.size(); ++t)
1161 : {
1162 : // get local solution at time point
1163 : LocalVector& locU = locTimeSeries.solution(t);
1164 0 : Eval.set_time_point(t);
1165 0 : number scaleStiff = vScaleStiff[t];
1166 :
1167 : // prepare element
1168 : try
1169 : {
1170 0 : Eval.prepare_elem(locU, elem, id, vCornerCoords, ind, false);
1171 : }
1172 0 : UG_CATCH_THROW("(instationary) AssembleLinear: Cannot prepare element.");
1173 :
1174 : // Assemble dM
1175 : try
1176 : {
1177 0 : tmpLocRhs = 0.0;
1178 0 : Eval.add_def_M_elem(tmpLocRhs, locU, elem, vCornerCoords, PT_INSTATIONARY);
1179 0 : locRhs.scale_append(-vScaleMass[t], tmpLocRhs);
1180 : }
1181 0 : UG_CATCH_THROW("(instationary) AssembleLinear: Cannot compute Jacobian (M).");
1182 :
1183 : // Assemble dA
1184 : try
1185 : {
1186 0 : if (scaleStiff != 0.0)
1187 : {
1188 0 : tmpLocRhs = 0.0;
1189 0 : Eval.add_def_A_elem(tmpLocRhs, locU, elem, vCornerCoords, PT_INSTATIONARY);
1190 0 : locRhs.scale_append(-scaleStiff, tmpLocRhs);
1191 : }
1192 : }
1193 0 : UG_CATCH_THROW("(instationary) AssembleLinear: Cannot compute Jacobian (A).");
1194 :
1195 : // Assemble rhs
1196 : try
1197 : {
1198 0 : if (scaleStiff != 0.0)
1199 : {
1200 0 : tmpLocRhs = 0.0;
1201 0 : Eval.add_rhs_elem(tmpLocRhs, elem, vCornerCoords, PT_INSTATIONARY);
1202 0 : locRhs.scale_append(scaleStiff, tmpLocRhs);
1203 : }
1204 : }
1205 0 : UG_CATCH_THROW("(instationary) AssembleLinear: Cannot compute Rhs.");
1206 : }
1207 :
1208 : // send local to global matrix & rhs
1209 : try{
1210 0 : if (!spAssTuner->matrix_is_const())
1211 0 : spAssTuner->add_local_mat_to_global(A, locA, dd);
1212 0 : spAssTuner->add_local_vec_to_global(rhs, locRhs, dd);
1213 : }
1214 0 : UG_CATCH_THROW("(instationary) AssembleLinear: Cannot add local vector/matrix.");
1215 : }
1216 :
1217 : // finish element loop
1218 : try
1219 : {
1220 0 : Eval.finish_elem_loop();
1221 : }
1222 0 : UG_CATCH_THROW("(instationary) AssembleLinear: Cannot finish element loop.");
1223 :
1224 0 : }
1225 0 : UG_CATCH_THROW("(instationary) AssembleLinear: Cannot create Data Evaluator.");
1226 : }
1227 :
1228 : ////////////////////////////////////////////////////////////////////////////////
1229 : // Assemble Rhs (of a stationary problem)
1230 : ////////////////////////////////////////////////////////////////////////////////
1231 :
1232 : public:
1233 : /**
1234 : * This function adds the contributions of all passed element discretizations
1235 : * on one given subset to the global Right-Hand Side. (This version processes
1236 : * elements in a given interval.)
1237 : *
1238 : * \param[in] vElemDisc element discretizations
1239 : * \param[in] spDomain domain
1240 : * \param[in] dd DoF Distribution
1241 : * \param[in] iterBegin element iterator
1242 : * \param[in] iterEnd element iterator
1243 : * \param[in] si subset index
1244 : * \param[in] bNonRegularGrid flag to indicate if non regular grid is used
1245 : * \param[in,out] rhs Right-hand side
1246 : * \param[in] u solution
1247 : * \param[in] spAssTuner assemble adapter
1248 : */
1249 : template <typename TElem, typename TIterator>
1250 : static void
1251 0 : AssembleRhs( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
1252 : ConstSmartPtr<domain_type> spDomain,
1253 : ConstSmartPtr<DoFDistribution> dd,
1254 : TIterator iterBegin,
1255 : TIterator iterEnd,
1256 : int si, bool bNonRegularGrid,
1257 : vector_type& rhs,
1258 : const vector_type& u,
1259 : ConstSmartPtr<AssemblingTuner<TAlgebra> > spAssTuner)
1260 : {
1261 : // check if there are any elements at all, otherwise return immediately
1262 0 : if(iterBegin == iterEnd) return;
1263 :
1264 : // reference object id
1265 : static const ReferenceObjectID id = geometry_traits<TElem>::REFERENCE_OBJECT_ID;
1266 :
1267 : // storage for corner coordinates
1268 0 : MathVector<domain_type::dim> vCornerCoords[TElem::NUM_VERTICES];
1269 :
1270 : // prepare for given elem discs
1271 : try
1272 : {
1273 0 : DataEvaluator<domain_type> Eval(RHS,
1274 : vElemDisc, dd->function_pattern(), bNonRegularGrid);
1275 :
1276 : // prepare loop
1277 0 : 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 0 : for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
1284 : {
1285 : // get Element
1286 0 : TElem* elem = *iter;
1287 :
1288 : // get corner coordinates
1289 0 : FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
1290 :
1291 : // check if elem is skipped from assembling
1292 0 : if(!spAssTuner->element_used(elem)) continue;
1293 :
1294 : // get global indices
1295 0 : dd->indices(elem, ind, Eval.use_hanging());
1296 :
1297 : // adapt local algebra
1298 0 : locU.resize(ind); locRhs.resize(ind);
1299 :
1300 : // read local values of u
1301 0 : GetLocalVector(locU, u);
1302 :
1303 : // prepare element
1304 : try
1305 : {
1306 0 : Eval.prepare_elem(locU, elem, id, vCornerCoords, ind);
1307 : }
1308 0 : UG_CATCH_THROW("AssembleRhs: Cannot prepare element.");
1309 :
1310 : // reset local algebra
1311 0 : locRhs = 0.0;
1312 :
1313 : // Assemble rhs
1314 : try
1315 : {
1316 0 : Eval.add_rhs_elem(locRhs, elem, vCornerCoords);
1317 : }
1318 0 : UG_CATCH_THROW("AssembleRhs: Cannot compute Rhs.");
1319 :
1320 : // send local to global rhs
1321 : try{
1322 0 : spAssTuner->add_local_vec_to_global(rhs, locRhs, dd);
1323 : }
1324 0 : UG_CATCH_THROW("AssembleRhs: Cannot add local vector.");
1325 : }
1326 :
1327 : // finish element loop
1328 : try
1329 : {
1330 0 : Eval.finish_elem_loop();
1331 : }
1332 0 : UG_CATCH_THROW("AssembleRhs: Cannot finish element loop.");
1333 :
1334 : }
1335 0 : UG_CATCH_THROW("AssembleRhs: Cannot create Data Evaluator.");
1336 : }
1337 :
1338 : ////////////////////////////////////////////////////////////////////////////////
1339 : // Assemble (instationary) Rhs
1340 : ////////////////////////////////////////////////////////////////////////////////
1341 :
1342 : public:
1343 : /**
1344 : * This function adds the contributions of all passed element discretizations
1345 : * on one given subset to the global Right-Hand Side in the time-dependent case.
1346 : * (This version processes elements in a given interval.)
1347 : *
1348 : * \param[in] vElemDisc element discretizations
1349 : * \param[in] spDomain domain
1350 : * \param[in] dd DoF Distribution
1351 : * \param[in] iterBegin element iterator
1352 : * \param[in] iterEnd element iterator
1353 : * \param[in] si subset index
1354 : * \param[in] bNonRegularGrid flag to indicate if non regular grid is used
1355 : * \param[in,out] rhs Right-hand side
1356 : * \param[in] vSol current and previous solutions
1357 : * \param[in] vScaleMass scaling factors for mass part
1358 : * \param[in] vScaleStiff scaling factors for stiffness part
1359 : * \param[in] spAssTuner assemble adapter
1360 : */
1361 : template <typename TElem, typename TIterator>
1362 : static void
1363 0 : AssembleRhs( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
1364 : ConstSmartPtr<domain_type> spDomain,
1365 : ConstSmartPtr<DoFDistribution> dd,
1366 : TIterator iterBegin,
1367 : TIterator iterEnd,
1368 : int si, bool bNonRegularGrid,
1369 : vector_type& rhs,
1370 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
1371 : const std::vector<number>& vScaleMass,
1372 : const std::vector<number>& vScaleStiff,
1373 : ConstSmartPtr<AssemblingTuner<TAlgebra> > spAssTuner)
1374 : {
1375 : // check if there are any elements at all, otherwise return immediately
1376 0 : if(iterBegin == iterEnd) return;
1377 :
1378 : // reference object id
1379 : static const ReferenceObjectID id = geometry_traits<TElem>::REFERENCE_OBJECT_ID;
1380 :
1381 : // storage for corner coordinates
1382 0 : MathVector<domain_type::dim> vCornerCoords[TElem::NUM_VERTICES];
1383 :
1384 : // check time scheme
1385 0 : if(vScaleMass.size() != vScaleStiff.size())
1386 0 : UG_THROW("(instationary) AssembleRhs: s_a and s_m must have same size.");
1387 :
1388 0 : if(vSol->size() < vScaleStiff.size())
1389 0 : 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 0 : locTimeSeries.read_times(vSol);
1396 :
1397 : // prepare for given elem discs
1398 : try
1399 : {
1400 0 : DataEvaluator<domain_type> Eval(MASS | STIFF | RHS,
1401 : vElemDisc, dd->function_pattern(), bNonRegularGrid,
1402 : &locTimeSeries, &vScaleMass, &vScaleStiff);
1403 :
1404 : // prepare loop
1405 0 : Eval.prepare_elem_loop(id, si);
1406 :
1407 : // local algebra
1408 : LocalIndices ind; LocalVector locRhs, tmpLocRhs;
1409 :
1410 : // Loop over all elements
1411 0 : for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
1412 : {
1413 : // get Element
1414 0 : TElem* elem = *iter;
1415 :
1416 : // get corner coordinates
1417 0 : FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
1418 :
1419 : // check if elem is skipped from assembling
1420 0 : if(!spAssTuner->element_used(elem)) continue;
1421 :
1422 : // get global indices
1423 0 : dd->indices(elem, ind, Eval.use_hanging());
1424 :
1425 : // adapt local algebra
1426 0 : locRhs.resize(ind); tmpLocRhs.resize(ind);
1427 :
1428 : // read local values of time series
1429 0 : locTimeSeries.read_values(vSol, ind);
1430 0 : Eval.set_time_point(0);
1431 :
1432 : // reset element contribution
1433 0 : locRhs = 0.0;
1434 :
1435 : /////////////////////
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 0 : Eval.prepare_elem(locU, elem, id, vCornerCoords, ind, false);
1445 : }
1446 0 : UG_CATCH_THROW("(instationary) AssembleRhs: Cannot prepare element.");
1447 :
1448 : // Assemble rhs
1449 : try
1450 : {
1451 0 : tmpLocRhs = 0.0;
1452 0 : Eval.add_rhs_elem(tmpLocRhs, elem, vCornerCoords, PT_INSTATIONARY);
1453 0 : locRhs.scale_append(vScaleStiff[0], tmpLocRhs);
1454 :
1455 0 : Eval.add_rhs_elem(locRhs, elem, vCornerCoords, PT_STATIONARY);
1456 : }
1457 0 : UG_CATCH_THROW("(instationary) AssembleRhs: Cannot compute Rhs.");
1458 :
1459 : ///////////////////
1460 : // old time steps
1461 :
1462 : // loop all old time points
1463 0 : for(size_t t = 1; t < vScaleStiff.size(); ++t)
1464 : {
1465 : // get local solution at time point
1466 : LocalVector& locU = locTimeSeries.solution(t);
1467 0 : Eval.set_time_point(t);
1468 :
1469 : // prepare element
1470 : try
1471 : {
1472 0 : Eval.prepare_elem(locU, elem, id, vCornerCoords, ind, false);
1473 : }
1474 0 : UG_CATCH_THROW("(instationary) AssembleRhs: Cannot prepare element.");
1475 :
1476 : // Assemble dM
1477 : try
1478 : {
1479 0 : tmpLocRhs = 0.0;
1480 0 : Eval.add_def_M_elem(tmpLocRhs, locU, elem, vCornerCoords, PT_INSTATIONARY);
1481 0 : locRhs.scale_append(-vScaleMass[t], tmpLocRhs);
1482 : }
1483 0 : UG_CATCH_THROW("(instationary) AssembleRhs: Cannot compute Jacobian (M).");
1484 :
1485 : // Assemble dA
1486 : try
1487 : {
1488 0 : tmpLocRhs = 0.0;
1489 0 : Eval.add_def_A_elem(tmpLocRhs, locU, elem, vCornerCoords, PT_INSTATIONARY);
1490 0 : locRhs.scale_append(-vScaleStiff[t], tmpLocRhs);
1491 : }
1492 0 : UG_CATCH_THROW("(instationary) AssembleRhs: Cannot compute Jacobian (A).");
1493 :
1494 : // Assemble rhs
1495 : try
1496 : {
1497 0 : tmpLocRhs = 0.0;
1498 0 : Eval.add_rhs_elem(tmpLocRhs, elem, vCornerCoords, PT_INSTATIONARY);
1499 0 : locRhs.scale_append(vScaleStiff[t], tmpLocRhs);
1500 : }
1501 0 : UG_CATCH_THROW("(instationary) AssembleRhs: Cannot compute Rhs.");
1502 : }
1503 :
1504 : // send local to global rhs
1505 : try{
1506 0 : spAssTuner->add_local_vec_to_global(rhs, locRhs, dd);
1507 : }
1508 0 : UG_CATCH_THROW("(instationary) AssembleRhs: Cannot add local vector.");
1509 : }
1510 :
1511 : // finish element loop
1512 : try
1513 : {
1514 0 : Eval.finish_elem_loop();
1515 : }
1516 0 : UG_CATCH_THROW("(instationary) AssembleRhs: Cannot finish element loop.");
1517 :
1518 : }
1519 0 : UG_CATCH_THROW("(instationary) AssembleRhs: Cannot create Data Evaluator.");
1520 : }
1521 :
1522 : // //////////////////////////////////////////////////////////////////////////////
1523 : // Prepare Timestep (for instationary problems)
1524 : // /////////////////////////////////////////////////////////////////////////////
1525 : public:
1526 : /**
1527 : * This function prepares the global discretization for a time-stepping scheme
1528 : * by calling the "prepare_timestep" methods of all passed element
1529 : * discretizations.
1530 : *
1531 : * \param[in] vElemDisc element discretizations
1532 : * \param[in] dd DoF Distribution
1533 : * \param[in] bNonRegularGrid flag to indicate if non regular grid is used
1534 : * \param[in] vSol current and previous solutions
1535 : * \param[in] spAssTuner assemble adapter
1536 : */
1537 : static void
1538 0 : PrepareTimestep
1539 : (
1540 : const std::vector<IElemDisc<domain_type>*>& vElemDisc,
1541 : ConstSmartPtr<DoFDistribution> dd,
1542 : bool bNonRegularGrid,
1543 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
1544 : number future_time,
1545 : ConstSmartPtr<AssemblingTuner<TAlgebra> > spAssTuner
1546 : )
1547 : {
1548 : // get current time and vector
1549 0 : const vector_type& u = *vSol->solution(0);
1550 :
1551 : // create local time series
1552 : LocalVectorTimeSeries locTimeSeries;
1553 0 : locTimeSeries.read_times(vSol);
1554 :
1555 : try
1556 : {
1557 0 : DataEvaluator<domain_type> Eval(NONE,
1558 : vElemDisc, dd->function_pattern(), bNonRegularGrid,
1559 : &locTimeSeries);
1560 0 : Eval.set_time_point(0);
1561 :
1562 : // prepare timestep
1563 : try
1564 : {
1565 : VectorProxy<vector_type> vp(u);
1566 0 : size_t algebra_index = bridge::AlgebraTypeIDProvider::instance().id<algebra_type>();
1567 0 : Eval.prepare_timestep(future_time, vSol->time(0), &vp, algebra_index);
1568 : }
1569 0 : UG_CATCH_THROW("(instationary) PrepareTimestep: Cannot prepare timestep.");
1570 :
1571 : }
1572 0 : UG_CATCH_THROW("(instationary) PrepareTimestep: Cannot create Data Evaluator.");
1573 0 : }
1574 :
1575 : ////////////////////////////////////////////////////////////////////////////////
1576 : // Prepare Timestep Elem (for instationary problems)
1577 : ////////////////////////////////////////////////////////////////////////////////
1578 : public:
1579 : /**
1580 : * This function prepares the global discretization for a time-stepping scheme
1581 : * by calling the "prepare_timestep_elem" methods of all passed element
1582 : * discretizations on one given subset.
1583 : * (This version processes elements in a given interval.)
1584 : *
1585 : * \param[in] vElemDisc element discretizations
1586 : * \param[in] spDomain domain
1587 : * \param[in] dd DoF Distribution
1588 : * \param[in] iterBegin element iterator
1589 : * \param[in] iterEnd element iterator
1590 : * \param[in] si subset index
1591 : * \param[in] bNonRegularGrid flag to indicate if non regular grid is used
1592 : * \param[in] vSol current and previous solutions
1593 : * \param[in] spAssTuner assemble adapter
1594 : */
1595 : template <typename TElem, typename TIterator>
1596 : static void
1597 0 : PrepareTimestepElem(const std::vector<IElemDisc<domain_type>*>& vElemDisc,
1598 : ConstSmartPtr<domain_type> spDomain,
1599 : ConstSmartPtr<DoFDistribution> dd,
1600 : TIterator iterBegin,
1601 : TIterator iterEnd,
1602 : int si, bool bNonRegularGrid,
1603 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
1604 : ConstSmartPtr<AssemblingTuner<TAlgebra> > spAssTuner)
1605 : {
1606 : // check if there are any elements at all, otherwise return immediately
1607 0 : if(iterBegin == iterEnd) return;
1608 :
1609 : // reference object id
1610 : static const ReferenceObjectID id = geometry_traits<TElem>::REFERENCE_OBJECT_ID;
1611 :
1612 : // storage for corner coordinates
1613 0 : MathVector<domain_type::dim> vCornerCoords[TElem::NUM_VERTICES];
1614 :
1615 : // get current time and vector
1616 0 : const vector_type& u = *vSol->solution(0);
1617 :
1618 : // create local time series
1619 : LocalVectorTimeSeries locTimeSeries;
1620 0 : locTimeSeries.read_times(vSol);
1621 :
1622 : try
1623 : {
1624 0 : DataEvaluator<domain_type> Eval(NONE,
1625 : vElemDisc, dd->function_pattern(), bNonRegularGrid,
1626 : &locTimeSeries);
1627 0 : Eval.set_time_point(0);
1628 :
1629 : // prepare element loop
1630 0 : Eval.prepare_elem_loop(id, si);
1631 :
1632 : // local algebra
1633 : LocalIndices ind; LocalVector locU;
1634 :
1635 : // Loop over all elements
1636 0 : for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
1637 : {
1638 : // get Element
1639 0 : TElem* elem = *iter;
1640 :
1641 : // get corner coordinates
1642 0 : FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
1643 :
1644 : // check if elem is skipped from assembling
1645 0 : if(!spAssTuner->element_used(elem)) continue;
1646 :
1647 : // get global indices
1648 0 : dd->indices(elem, ind, Eval.use_hanging());
1649 :
1650 : // adapt local algebra
1651 0 : locU.resize(ind);
1652 :
1653 : // read local values of u
1654 0 : GetLocalVector(locU, u);
1655 :
1656 : // read local values of time series
1657 0 : if(Eval.time_series_needed())
1658 0 : locTimeSeries.read_values(vSol, ind);
1659 :
1660 : // prepare timestep
1661 : try
1662 : {
1663 0 : Eval.prepare_timestep_elem(vSol->time(0), locU, elem, vCornerCoords);
1664 : }
1665 0 : UG_CATCH_THROW("(instationary) PrepareTimestep: Cannot prepare timestep.");
1666 : }
1667 :
1668 : // finish element loop
1669 : try
1670 : {
1671 0 : Eval.finish_elem_loop();
1672 : }
1673 0 : UG_CATCH_THROW("AssembleMassMatrix: Cannot finish element loop.");
1674 :
1675 : }
1676 0 : UG_CATCH_THROW("(instationary) PrepareTimestep: Cannot create Data Evaluator.");
1677 : }
1678 :
1679 : // //////////////////////////////////////////////////////////////////////////////
1680 : // Finish Timestep (for instationary problems)
1681 : // /////////////////////////////////////////////////////////////////////////////
1682 : public:
1683 : /**
1684 : * This function finishes the global discretization for a time-stepping scheme
1685 : * by calling the "finish_timestep" methods of all passed element
1686 : * discretizations.
1687 : *
1688 : * \param[in] vElemDisc element discretizations
1689 : * \param[in] dd DoF Distribution
1690 : * \param[in] bNonRegularGrid flag to indicate if non regular grid is used
1691 : * \param[in] vSol current and previous solutions
1692 : * \param[in] spAssTuner assemble adapter
1693 : */
1694 : static void
1695 0 : FinishTimestep
1696 : (
1697 : const std::vector<IElemDisc<domain_type>*>& vElemDisc,
1698 : ConstSmartPtr<DoFDistribution> dd,
1699 : bool bNonRegularGrid,
1700 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
1701 : ConstSmartPtr<AssemblingTuner<TAlgebra> > spAssTuner
1702 : )
1703 : {
1704 : // get current time and vector
1705 0 : const vector_type& u = *vSol->solution(0);
1706 :
1707 : // create local time series
1708 : LocalVectorTimeSeries locTimeSeries;
1709 0 : locTimeSeries.read_times(vSol);
1710 :
1711 : try
1712 : {
1713 0 : DataEvaluator<domain_type> Eval(NONE,
1714 : vElemDisc, dd->function_pattern(), bNonRegularGrid,
1715 : &locTimeSeries);
1716 0 : Eval.set_time_point(0);
1717 :
1718 : // finish time step
1719 : try
1720 : {
1721 : VectorProxy<vector_type> vp(u);
1722 0 : size_t algebra_index = bridge::AlgebraTypeIDProvider::instance().id<algebra_type>();
1723 0 : Eval.finish_timestep(vSol->time(0), &vp, algebra_index);
1724 : }
1725 0 : UG_CATCH_THROW("(instationary) FinishTimestep: Cannot prepare time step.");
1726 :
1727 : }
1728 0 : UG_CATCH_THROW("(instationary) FinishTimestep: Cannot create Data Evaluator.");
1729 0 : }
1730 :
1731 : ////////////////////////////////////////////////////////////////////////////////
1732 : // Finish Timestep Elem (for instationary problems)
1733 : ////////////////////////////////////////////////////////////////////////////////
1734 :
1735 : public:
1736 : /**
1737 : * This function finalizes the global discretization in a time-stepping scheme
1738 : * by calling the "finish_timestep_elem" methods of all passed element
1739 : * discretizations on one given subset.
1740 : * (This version processes elements in a given interval.)
1741 : *
1742 : * \param[in] vElemDisc element discretizations
1743 : * \param[in] dd DoF Distribution
1744 : * \param[in] iterBegin element iterator
1745 : * \param[in] iterEnd element iterator
1746 : * \param[in] si subset index
1747 : * \param[in] bNonRegularGrid flag to indicate if non regular grid is used
1748 : * \param[in] vSol current and previous solutions
1749 : * \param[in] spAssTuner assemble adapter
1750 : */
1751 : template <typename TElem, typename TIterator>
1752 : static void
1753 0 : FinishTimestepElem(const std::vector<IElemDisc<domain_type>*>& vElemDisc,
1754 : ConstSmartPtr<domain_type> spDomain,
1755 : ConstSmartPtr<DoFDistribution> dd,
1756 : TIterator iterBegin,
1757 : TIterator iterEnd,
1758 : int si, bool bNonRegularGrid,
1759 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
1760 : ConstSmartPtr<AssemblingTuner<TAlgebra> > spAssTuner)
1761 : {
1762 : // check if there are any elements at all, otherwise return immediately
1763 0 : if(iterBegin == iterEnd) return;
1764 :
1765 : // reference object id
1766 : static const ReferenceObjectID id = geometry_traits<TElem>::REFERENCE_OBJECT_ID;
1767 :
1768 : // storage for corner coordinates
1769 0 : MathVector<domain_type::dim> vCornerCoords[TElem::NUM_VERTICES];
1770 :
1771 : // get current time and vector
1772 0 : const vector_type& u = *vSol->solution(0);
1773 :
1774 : // create local time series
1775 : LocalVectorTimeSeries locTimeSeries;
1776 0 : locTimeSeries.read_times(vSol);
1777 :
1778 : // prepare for given elem discs
1779 : try
1780 : {
1781 0 : DataEvaluator<domain_type> Eval(NONE,
1782 : vElemDisc, dd->function_pattern(), bNonRegularGrid,
1783 : &locTimeSeries);
1784 0 : Eval.set_time_point(0);
1785 :
1786 : // prepare loop
1787 0 : Eval.prepare_elem_loop(id, si);
1788 :
1789 : // local algebra
1790 : LocalIndices ind; LocalVector locU;
1791 :
1792 : // Loop over all elements
1793 0 : for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
1794 : {
1795 : // get Element
1796 0 : TElem* elem = *iter;
1797 :
1798 : // get corner coordinates
1799 0 : FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
1800 :
1801 : // check if elem is skipped from assembling
1802 0 : if(!spAssTuner->element_used(elem)) continue;
1803 :
1804 : // get global indices
1805 0 : dd->indices(elem, ind, Eval.use_hanging());
1806 :
1807 : // adapt local algebra
1808 0 : locU.resize(ind);
1809 :
1810 : // read local values of u
1811 0 : GetLocalVector(locU, u);
1812 :
1813 : // read local values of time series
1814 0 : if(Eval.time_series_needed())
1815 0 : locTimeSeries.read_values(vSol, ind);
1816 :
1817 : // finish timestep
1818 : try
1819 : {
1820 0 : Eval.finish_timestep_elem(locTimeSeries.time(0), locU, elem, vCornerCoords);
1821 : }
1822 0 : UG_CATCH_THROW("(instationary) FinishTimestepElem: Cannot finish timestep.");
1823 : }
1824 :
1825 : // finish element loop
1826 : try
1827 : {
1828 0 : Eval.finish_elem_loop();
1829 : }
1830 0 : UG_CATCH_THROW("AssembleMassMatrix: Cannot finish element loop.");
1831 :
1832 : }
1833 0 : UG_CATCH_THROW("(instationary) FinishTimestepElem: Cannot create Data Evaluator");
1834 : }
1835 :
1836 : ////////////////////////////////////////////////////////////////////////////////
1837 : // Init. all exports (an optional operation, to use the exports for plotting etc.)
1838 : ////////////////////////////////////////////////////////////////////////////////
1839 :
1840 : public:
1841 : /**
1842 : * This function finalizes the global discretization in a time-stepping scheme
1843 : * by calling the "finish_timestep_elem" methods of all passed element
1844 : * discretizations on one given subset.
1845 : * (This version processes elements in a given interval.)
1846 : *
1847 : * \param[in] vElemDisc element discretizations
1848 : * \param[in] dd DoF Distribution
1849 : * \param[in] iterBegin element iterator
1850 : * \param[in] iterEnd element iterator
1851 : * \param[in] si subset index
1852 : * \param[in] bNonRegularGrid flag to indicate if non regular grid is used
1853 : * \param[in] bAsTimeDependent flag to simulate the instationary case for the initialization
1854 : */
1855 : template <typename TElem, typename TIterator>
1856 : static void
1857 0 : InitAllExports(const std::vector<IElemDisc<domain_type>*>& vElemDisc,
1858 : ConstSmartPtr<DoFDistribution> dd,
1859 : TIterator iterBegin,
1860 : TIterator iterEnd,
1861 : int si, bool bNonRegularGrid, bool bAsTimeDependent)
1862 : {
1863 : // check if there are any elements at all, otherwise return immediately
1864 0 : if(iterBegin == iterEnd) return;
1865 :
1866 : // reference object id
1867 : static const ReferenceObjectID id = geometry_traits<TElem>::REFERENCE_OBJECT_ID;
1868 :
1869 : // dummy local time series (only to simulate the instationary case for the initialization)
1870 : LocalVectorTimeSeries locTimeSeries;
1871 :
1872 : // prepare for given elem discs
1873 : try
1874 : {
1875 0 : DataEvaluator<domain_type> Eval(MASS | STIFF | RHS,
1876 : vElemDisc, dd->function_pattern(), bNonRegularGrid,
1877 : bAsTimeDependent? &locTimeSeries : NULL);
1878 0 : Eval.set_time_point(0);
1879 :
1880 : // prepare loop
1881 0 : Eval.prepare_elem_loop(id, si);
1882 :
1883 :
1884 : // finish element loop
1885 0 : Eval.finish_elem_loop();
1886 :
1887 : }
1888 0 : UG_CATCH_THROW("InitAllExports: Data Evaluator failed");
1889 : }
1890 :
1891 : ////////////////////////////////////////////////////////////////////////////////
1892 : // Error estimator (stationary)
1893 : ////////////////////////////////////////////////////////////////////////////////
1894 :
1895 : public:
1896 : /**
1897 : * This function assembles the error estimator associated with all the
1898 : * element discretizations in the internal data structure for one given subset.
1899 : * (This version processes elements in a given interval.)
1900 : *
1901 : * \param[in] vElemDisc element discretizations
1902 : * \param[in] spDomain domain
1903 : * \param[in] dd DoF Distribution
1904 : * \param[in] iterBegin element iterator
1905 : * \param[in] iterEnd element iterator
1906 : * \param[in] si subset index
1907 : * \param[in] bNonRegularGrid flag to indicate if non regular grid is used
1908 : * \param[in] u solution
1909 : */
1910 : template <typename TElem, typename TIterator>
1911 : static void
1912 0 : AssembleErrorEstimator
1913 : (
1914 : const std::vector<IElemError<domain_type>*>& vElemDisc,
1915 : ConstSmartPtr<domain_type> spDomain,
1916 : ConstSmartPtr<DoFDistribution> dd,
1917 : TIterator iterBegin,
1918 : TIterator iterEnd,
1919 : int si,
1920 : bool bNonRegularGrid,
1921 : const vector_type& u
1922 : )
1923 : {
1924 : // check if there are any elements at all, otherwise return immediately
1925 0 : if(iterBegin == iterEnd) return;
1926 :
1927 : // reference object id
1928 : static const ReferenceObjectID id = geometry_traits<TElem>::REFERENCE_OBJECT_ID;
1929 :
1930 : // storage for corner coordinates
1931 0 : MathVector<domain_type::dim> vCornerCoords[TElem::NUM_VERTICES];
1932 :
1933 : // prepare for given elem discs
1934 : try
1935 : {
1936 0 : ErrorEvaluator<domain_type> Eval(STIFF | RHS,
1937 : vElemDisc, dd->function_pattern(), bNonRegularGrid);
1938 :
1939 : // prepare element loop
1940 0 : Eval.prepare_err_est_elem_loop(id, si);
1941 :
1942 : // local indices and local algebra
1943 : LocalIndices ind; LocalVector locU;
1944 :
1945 : // Loop over all elements
1946 0 : for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
1947 : {
1948 : // get Element
1949 : TElem* elem = *iter;
1950 :
1951 : // get corner coordinates
1952 0 : FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
1953 :
1954 : // get global indices
1955 0 : dd->indices(elem, ind, Eval.use_hanging());
1956 :
1957 : // adapt local algebra
1958 0 : locU.resize(ind);
1959 :
1960 : // read local values of u
1961 0 : GetLocalVector(locU, u);
1962 :
1963 : // prepare element
1964 : try
1965 : {
1966 0 : Eval.prepare_err_est_elem(locU, elem, vCornerCoords, ind, false);
1967 : }
1968 0 : UG_CATCH_THROW("(stationary) AssembleRhs: Cannot prepare element.");
1969 :
1970 : // assemble stiffness part
1971 : try
1972 : {
1973 0 : Eval.compute_err_est_A_elem(locU, elem, vCornerCoords, ind);
1974 : }
1975 0 : UG_CATCH_THROW("AssembleErrorEstimator: Cannot assemble the error estimator for stiffness part.");
1976 :
1977 : // assemble rhs part
1978 : try
1979 : {
1980 0 : Eval.compute_err_est_rhs_elem(elem, vCornerCoords, ind);
1981 : }
1982 0 : UG_CATCH_THROW("AssembleErrorEstimator: Cannot assemble the error estimator for stiffness part.");
1983 :
1984 : }
1985 :
1986 : // finish element loop
1987 : try
1988 : {
1989 0 : Eval.finish_err_est_elem_loop();
1990 : }
1991 0 : UG_CATCH_THROW("AssembleErrorEstimator: Cannot finish element loop.");
1992 :
1993 : }
1994 0 : UG_CATCH_THROW("AssembleErrorEstimator: Cannot create Data Evaluator.");
1995 : }
1996 :
1997 : ////////////////////////////////////////////////////////////////////////////////
1998 : // Error estimator (for time-dependent problems)
1999 : ////////////////////////////////////////////////////////////////////////////////
2000 :
2001 : public:
2002 : /**
2003 : * This function assembles the error estimator associated with all the
2004 : * element discretizations in the internal data structure for one given subset.
2005 : * (This version processes elements in a given interval.)
2006 : *
2007 : * \param[in] vElemDisc element discretizations
2008 : * \param[in] spDomain domain
2009 : * \param[in] dd DoF Distribution
2010 : * \param[in] iterBegin element iterator
2011 : * \param[in] iterEnd element iterator
2012 : * \param[in] si subset index
2013 : * \param[in] bNonRegularGrid flag to indicate if non regular grid is used
2014 : * \param[in] vScaleMass scaling factors for mass part
2015 : * \param[in] vScaleStiff scaling factors for stiffness part
2016 : * \param[in] vSol solution
2017 : */
2018 : template <typename TElem, typename TIterator>
2019 : static void
2020 0 : AssembleErrorEstimator
2021 : (
2022 : const std::vector<IElemError<domain_type>*>& vElemDisc,
2023 : ConstSmartPtr<domain_type> spDomain,
2024 : ConstSmartPtr<DoFDistribution> dd,
2025 : TIterator iterBegin,
2026 : TIterator iterEnd,
2027 : int si,
2028 : bool bNonRegularGrid,
2029 : const std::vector<number>& vScaleMass,
2030 : const std::vector<number>& vScaleStiff,
2031 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol
2032 : )
2033 : {
2034 : // check if there are any elements at all, otherwise return immediately
2035 0 : if (iterBegin == iterEnd) return;
2036 :
2037 : // reference object id
2038 : static const ReferenceObjectID id = geometry_traits<TElem>::REFERENCE_OBJECT_ID;
2039 :
2040 : // storage for corner coordinates
2041 0 : MathVector<domain_type::dim> vCornerCoords[TElem::NUM_VERTICES];
2042 :
2043 : // check time scheme
2044 0 : if(vScaleMass.size() != vScaleStiff.size())
2045 0 : UG_THROW("(instationary) AssembleErrorEstimator: s_a and s_m must have same size.");
2046 :
2047 0 : if(vSol->size() < vScaleStiff.size())
2048 0 : UG_THROW("(instationary) AssembleErrorEstimator: Time stepping scheme needs at "
2049 : "least "<<vScaleStiff.size()<<" time steps, but only "<<
2050 : vSol->size() << " passed.");
2051 :
2052 : // create local time series
2053 : LocalVectorTimeSeries locTimeSeries;
2054 0 : locTimeSeries.read_times(vSol);
2055 :
2056 : // prepare for given elem discs
2057 : try
2058 : {
2059 0 : ErrorEvaluator<domain_type> Eval(MASS | STIFF | RHS,
2060 : vElemDisc, dd->function_pattern(), bNonRegularGrid,
2061 : &locTimeSeries, &vScaleMass, &vScaleStiff);
2062 :
2063 : // prepare element loop
2064 0 : Eval.prepare_err_est_elem_loop(id, si);
2065 :
2066 : // local indices and local algebra
2067 : LocalIndices ind;
2068 :
2069 : // loop over all elements
2070 0 : for (TIterator iter = iterBegin; iter != iterEnd; ++iter)
2071 : {
2072 : // get Element
2073 : TElem* elem = *iter;
2074 :
2075 : // get corner coordinates
2076 0 : FillCornerCoordinates(vCornerCoords, *elem, *spDomain);
2077 :
2078 : // get global indices
2079 0 : dd->indices(elem, ind, Eval.use_hanging());
2080 :
2081 : // read local values of time series
2082 0 : locTimeSeries.read_values(vSol, ind);
2083 :
2084 : // loop all time points and assemble them
2085 0 : for (std::size_t t = 0; t < vScaleStiff.size(); ++t)
2086 : {
2087 : // get local solution at timepoint
2088 : LocalVector& locU = locTimeSeries.solution(t);
2089 0 : Eval.set_time_point(t);
2090 :
2091 : // prepare element
2092 : try
2093 : {
2094 0 : Eval.prepare_err_est_elem(locU, elem, vCornerCoords, ind, false);
2095 : }
2096 0 : UG_CATCH_THROW("AssembleErrorEstimator: Cannot prepare element.");
2097 :
2098 : // assemble stiffness part
2099 : try
2100 : {
2101 0 : Eval.compute_err_est_A_elem(locU, elem, vCornerCoords, ind, vScaleMass[t], vScaleStiff[t]);
2102 : }
2103 0 : UG_CATCH_THROW("AssembleErrorEstimator: Cannot assemble the error estimator for stiffness part.");
2104 :
2105 : // assemble mass part
2106 : try
2107 : {
2108 0 : Eval.compute_err_est_M_elem(locU, elem, vCornerCoords, ind, vScaleMass[t], vScaleStiff[t]);
2109 : }
2110 0 : UG_CATCH_THROW("AssembleErrorEstimator: Cannot assemble the error estimator for stiffness part.");
2111 :
2112 : // assemble rhs part
2113 : try
2114 : {
2115 0 : Eval.compute_err_est_rhs_elem(elem, vCornerCoords, ind, vScaleMass[t], vScaleStiff[t]);
2116 : }
2117 0 : UG_CATCH_THROW("AssembleErrorEstimator: Cannot assemble the error estimator for stiffness part.");
2118 : }
2119 : }
2120 :
2121 : // finish element loop
2122 : try
2123 : {
2124 0 : Eval.finish_err_est_elem_loop();
2125 : }
2126 0 : UG_CATCH_THROW("AssembleErrorEstimator: Cannot finish element loop.");
2127 :
2128 : }
2129 0 : UG_CATCH_THROW("AssembleErrorEstimator: Cannot create Data Evaluator.");
2130 : }
2131 :
2132 : }; // class StdGlobAssembler
2133 :
2134 : } // end namespace ug
2135 :
2136 :
2137 : #endif /* __H__UG__LIB_DISC__SPATIAL_DISC__ELEM_DISC__ELEM_DISC_ASSEMBLE_UTIL__ */
|
