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__DOMAIN_DISC_IMPL__
34 : #define __H__UG__LIB_DISC__SPATIAL_DISC__DOMAIN_DISC_IMPL__
35 :
36 : #include "common/profiler/profiler.h"
37 : #include "domain_disc.h"
38 : #include "lib_disc/common/groups_util.h"
39 : #include "lib_disc/function_spaces/error_indicator_util.h"
40 : #include "lib_disc/spatial_disc/subset_assemble_util.h"
41 : #ifdef UG_PARALLEL
42 : #include "lib_disc/parallelization/parallelization_util.h"
43 : #endif
44 : #include "lib_grid/algorithms/debug_util.h"
45 :
46 : namespace ug{
47 :
48 : template <class TElemDisc>
49 : static void prep_assemble_loop(std::vector<TElemDisc*> vElemDisc)
50 : {
51 0 : for(size_t i = 0; i < vElemDisc.size(); ++i)
52 : {
53 0 : vElemDisc[i]->prep_assemble_loop();
54 : }
55 : }
56 :
57 : template <class TElemDisc>
58 : static void post_assemble_loop(std::vector<TElemDisc*> vElemDisc)
59 : {
60 0 : for(size_t i = 0; i < vElemDisc.size(); ++i)
61 : {
62 0 : vElemDisc[i]->post_assemble_loop();
63 : }
64 : }
65 :
66 :
67 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
68 0 : void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::update_elem_discs()
69 : {
70 : // check Approximation space
71 0 : if(!m_spApproxSpace.valid())
72 0 : UG_THROW("DomainDiscretization: Before using the "
73 : "DomainDiscretization an ApproximationSpace must be set to it. "
74 : "Please use DomainDiscretization:set_approximation_space to "
75 : "set an appropriate Space.");
76 :
77 : // set approximation space and extract IElemDiscs
78 : m_vElemDisc.clear();
79 0 : for(size_t i = 0; i < m_vDomainElemDisc.size(); ++i)
80 : {
81 0 : m_vDomainElemDisc[i]->set_approximation_space(m_spApproxSpace);
82 :
83 0 : if(!(m_spAssTuner->elem_disc_type_enabled(m_vDomainElemDisc[i]->type()))) continue;
84 0 : m_vElemDisc.push_back(m_vDomainElemDisc[i].get());
85 : }
86 0 : }
87 :
88 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
89 0 : void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::update_elem_errors()
90 : {
91 : // check Approximation space
92 0 : if(!m_spApproxSpace.valid())
93 0 : UG_THROW("DomainDiscretization: Before using the "
94 : "DomainDiscretization an ApproximationSpace must be set to it. "
95 : "Please use DomainDiscretization:set_approximation_space to "
96 : "set an appropriate Space.");
97 :
98 : // set approximation space and extract IElemDiscs
99 : m_vElemError.clear();
100 0 : for(size_t i = 0; i < m_vDomainElemError.size(); ++i)
101 : {
102 0 : m_vDomainElemError[i]->set_approximation_space(m_spApproxSpace);
103 :
104 0 : if(!(m_spAssTuner->elem_disc_type_enabled(m_vDomainElemError[i]->type()))) continue;
105 0 : m_vElemError.push_back(m_vDomainElemError[i].get());
106 : }
107 :
108 0 : }
109 :
110 :
111 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
112 0 : void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::update_constraints()
113 : {
114 : // check Approximation space
115 0 : if(!m_spApproxSpace.valid())
116 0 : UG_THROW("DomainDiscretization: Before using the "
117 : "DomainDiscretization an ApproximationSpace must be set to it. "
118 : "Please use DomainDiscretization:set_approximation_space to "
119 : "set an appropriate Space.");
120 :
121 :
122 0 : for(size_t i = 0; i < m_vConstraint.size(); ++i)
123 0 : m_vConstraint[i]->set_approximation_space(m_spApproxSpace);
124 0 : }
125 :
126 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
127 : void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::update_disc_items()
128 : {
129 0 : update_elem_discs();
130 0 : update_constraints();
131 0 : }
132 :
133 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
134 : void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::update_error_items()
135 : {
136 0 : update_elem_errors();
137 0 : update_constraints();
138 0 : }
139 :
140 : ///////////////////////////////////////////////////////////////////////////////
141 : // Mass Matrix
142 : ///////////////////////////////////////////////////////////////////////////////
143 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
144 0 : void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::
145 : assemble_mass_matrix(matrix_type& M, const vector_type& u,
146 : ConstSmartPtr<DoFDistribution> dd)
147 : {
148 : PROFILE_FUNC_GROUP("discretization");
149 : // update the elem discs
150 : update_disc_items();
151 0 : prep_assemble_loop(m_vElemDisc);
152 :
153 : // reset matrix to zero and resize
154 0 : m_spAssTuner->resize(dd, M);
155 :
156 : // Union of Subsets
157 0 : SubsetGroup unionSubsets;
158 : std::vector<SubsetGroup> vSSGrp;
159 :
160 : // create list of all subsets
161 : try{
162 0 : CreateSubsetGroups(vSSGrp, unionSubsets, m_vElemDisc, dd->subset_handler());
163 0 : }UG_CATCH_THROW("'DomainDiscretization': Can not create Subset Groups and Union.");
164 :
165 : // loop subsets
166 0 : for(size_t i = 0; i < unionSubsets.size(); ++i)
167 : {
168 : // get subset
169 : const int si = unionSubsets[i];
170 :
171 : // get dimension of the subset
172 0 : const int dim = DimensionOfSubset(*dd->subset_handler(), si);
173 :
174 : // request if subset is regular grid
175 0 : bool bNonRegularGrid = !unionSubsets.regular_grid(i);
176 :
177 : // overrule by regular grid if required
178 0 : if(m_spAssTuner->regular_grid_forced()) bNonRegularGrid = false;
179 :
180 : // Elem Disc on the subset
181 : std::vector<IElemDisc<TDomain>*> vSubsetElemDisc;
182 :
183 : // get all element discretizations that work on the subset
184 0 : GetElemDiscOnSubset(vSubsetElemDisc, m_vElemDisc, vSSGrp, si);
185 :
186 : // assemble on suitable elements
187 : try
188 : {
189 0 : switch(dim)
190 : {
191 0 : case 0:
192 : this->template AssembleMassMatrix<RegularVertex>
193 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, M, u);
194 0 : break;
195 0 : case 1:
196 : this->template AssembleMassMatrix<RegularEdge>
197 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, M, u);
198 : // When assembling over lower-dim manifolds that contain hanging nodes:
199 : this->template AssembleMassMatrix<ConstrainingEdge>
200 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, M, u);
201 0 : break;
202 0 : case 2:
203 : this->template AssembleMassMatrix<Triangle>
204 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, M, u);
205 : this->template AssembleMassMatrix<Quadrilateral>
206 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, M, u);
207 : // When assembling over lower-dim manifolds that contain hanging nodes:
208 : this->template AssembleMassMatrix<ConstrainingTriangle>
209 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, M, u);
210 : this->template AssembleMassMatrix<ConstrainingQuadrilateral>
211 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, M, u);
212 0 : break;
213 0 : case 3:
214 : this->template AssembleMassMatrix<Tetrahedron>
215 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, M, u);
216 : this->template AssembleMassMatrix<Pyramid>
217 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, M, u);
218 : this->template AssembleMassMatrix<Prism>
219 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, M, u);
220 : this->template AssembleMassMatrix<Hexahedron>
221 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, M, u);
222 : this->template AssembleMassMatrix<Octahedron>
223 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, M, u);
224 0 : break;
225 0 : default:
226 0 : UG_THROW("DomainDiscretization::assemble_mass_matrix:"
227 : "Dimension "<<dim<<" (subset="<<si<<") not supported.");
228 : }
229 : }
230 0 : UG_CATCH_THROW("DomainDiscretization::assemble_mass_matrix:"
231 : " Assembling of elements of Dimension " << dim << " in "
232 : " subset "<<si<< " failed.");
233 : }
234 :
235 : // post process
236 : try{
237 0 : for(int type = 1; type < CT_ALL; type = type << 1){
238 0 : if(!(m_spAssTuner->constraint_type_enabled(type))) continue;
239 0 : for(size_t i = 0; i < m_vConstraint.size(); ++i)
240 0 : if(m_vConstraint[i]->type() & type)
241 : {
242 0 : m_vConstraint[i]->set_ass_tuner(m_spAssTuner);
243 0 : m_vConstraint[i]->adjust_jacobian(M, u, dd, type);
244 : }
245 : }
246 0 : post_assemble_loop(m_vElemDisc);
247 0 : }UG_CATCH_THROW("DomainDiscretization::assemble_mass_matrix:"
248 : " Cannot execute post process.");
249 :
250 : // Remember parallel storage type
251 : #ifdef UG_PARALLEL
252 : M.set_storage_type(PST_ADDITIVE);
253 : M.set_layouts(dd->layouts());
254 : #endif
255 0 : }
256 :
257 : /**
258 : * This function adds the contributions of all passed element discretizations
259 : * on one given subset to the global Mass matrix.
260 : *
261 : * \param[in] vElemDisc element discretizations
262 : * \param[in] dd DoF Distribution
263 : * \param[in] si subset index
264 : * \param[in] bNonRegularGrid flag to indicate if non regular grid is used
265 : * \param[in,out] M Mass matrix
266 : * \param[in] u solution
267 : */
268 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
269 : template <typename TElem>
270 0 : void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::
271 : AssembleMassMatrix( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
272 : ConstSmartPtr<DoFDistribution> dd,
273 : int si, bool bNonRegularGrid,
274 : matrix_type& M,
275 : const vector_type& u)
276 : {
277 : // check if only some elements are selected
278 0 : if(m_spAssTuner->selected_elements_used())
279 : {
280 : std::vector<TElem*> vElem;
281 0 : m_spAssTuner->collect_selected_elements(vElem, dd, si);
282 :
283 : // assembling is carried out only over those elements
284 : // which are selected and in subset si
285 : gass_type::template AssembleMassMatrix<TElem>
286 0 : (vElemDisc, m_spApproxSpace->domain(), dd, vElem.begin(), vElem.end(), si,
287 : bNonRegularGrid, M, u, m_spAssTuner);
288 0 : }
289 : else
290 : {
291 : // general case: assembling over all elements in subset si
292 : gass_type::template AssembleMassMatrix<TElem>
293 0 : (vElemDisc, m_spApproxSpace->domain(), dd,
294 : dd->template begin<TElem>(si), dd->template end<TElem>(si), si,
295 : bNonRegularGrid, M, u, m_spAssTuner);
296 : }
297 0 : }
298 :
299 : ///////////////////////////////////////////////////////////////////////////////
300 : // Stiffness Matrix
301 : ///////////////////////////////////////////////////////////////////////////////
302 :
303 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
304 0 : void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::
305 : assemble_stiffness_matrix(matrix_type& A, const vector_type& u,
306 : ConstSmartPtr<DoFDistribution> dd)
307 : {
308 : PROFILE_FUNC_GROUP("discretization");
309 : // update the elem discs
310 : update_disc_items();
311 0 : prep_assemble_loop(m_vElemDisc);
312 :
313 : // reset matrix to zero and resize
314 0 : m_spAssTuner->resize(dd, A);
315 :
316 : // Union of Subsets
317 0 : SubsetGroup unionSubsets;
318 : std::vector<SubsetGroup> vSSGrp;
319 :
320 : // create list of all subsets
321 : try{
322 0 : CreateSubsetGroups(vSSGrp, unionSubsets, m_vElemDisc, dd->subset_handler());
323 0 : }UG_CATCH_THROW("'DomainDiscretization': Can not create Subset Groups and Union.");
324 :
325 : // loop subsets
326 0 : for(size_t i = 0; i < unionSubsets.size(); ++i)
327 : {
328 : // get subset
329 : const int si = unionSubsets[i];
330 :
331 : // get dimension of the subset
332 0 : const int dim = DimensionOfSubset(*dd->subset_handler(), si);
333 :
334 : // request if subset is regular grid
335 0 : bool bNonRegularGrid = !unionSubsets.regular_grid(i);
336 :
337 : // overrule by regular grid if required
338 0 : if(m_spAssTuner->regular_grid_forced()) bNonRegularGrid = false;
339 :
340 : // Elem Disc on the subset
341 : std::vector<IElemDisc<TDomain>*> vSubsetElemDisc;
342 :
343 : // get all element discretizations that work on the subset
344 0 : GetElemDiscOnSubset(vSubsetElemDisc, m_vElemDisc, vSSGrp, si);
345 :
346 : // assemble on suitable elements
347 : try
348 : {
349 0 : switch(dim)
350 : {
351 0 : case 0:
352 : this->template AssembleStiffnessMatrix<RegularVertex>
353 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, A, u);
354 0 : break;
355 0 : case 1:
356 : this->template AssembleStiffnessMatrix<RegularEdge>
357 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, A, u);
358 : // When assembling over lower-dim manifolds that contain hanging nodes:
359 : this->template AssembleStiffnessMatrix<ConstrainingEdge>
360 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, A, u);
361 0 : break;
362 0 : case 2:
363 : this->template AssembleStiffnessMatrix<Triangle>
364 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, A, u);
365 : this->template AssembleStiffnessMatrix<Quadrilateral>
366 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, A, u);
367 : // When assembling over lower-dim manifolds that contain hanging nodes:
368 : this->template AssembleStiffnessMatrix<ConstrainingTriangle>
369 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, A, u);
370 : this->template AssembleStiffnessMatrix<ConstrainingQuadrilateral>
371 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, A, u);
372 0 : break;
373 0 : case 3:
374 : this->template AssembleStiffnessMatrix<Tetrahedron>
375 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, A, u);
376 : this->template AssembleStiffnessMatrix<Pyramid>
377 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, A, u);
378 : this->template AssembleStiffnessMatrix<Prism>
379 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, A, u);
380 : this->template AssembleStiffnessMatrix<Hexahedron>
381 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, A, u);
382 : this->template AssembleStiffnessMatrix<Octahedron>
383 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, A, u);
384 0 : break;
385 0 : default:
386 0 : UG_THROW("DomainDiscretization::assemble_stiffness_matrix:"
387 : "Dimension "<<dim<<" (subset="<<si<<") not supported.");
388 : }
389 : }
390 0 : UG_CATCH_THROW("DomainDiscretization::assemble_stiffness_matrix:"
391 : " Assembling of elements of Dimension " << dim << " in "
392 : " subset "<<si<< " failed.");
393 : }
394 :
395 : // post process
396 : try{
397 0 : for(int type = 1; type < CT_ALL; type = type << 1){
398 0 : if(!(m_spAssTuner->constraint_type_enabled(type))) continue;
399 0 : for(size_t i = 0; i < m_vConstraint.size(); ++i)
400 0 : if(m_vConstraint[i]->type() & type)
401 : {
402 0 : m_vConstraint[i]->set_ass_tuner(m_spAssTuner);
403 0 : m_vConstraint[i]->adjust_jacobian(A, u, dd, type);
404 : }
405 : }
406 0 : post_assemble_loop(m_vElemDisc);
407 0 : }UG_CATCH_THROW("DomainDiscretization::assemble_stiffness_matrix:"
408 : " Cannot execute post process.");
409 :
410 : // Remember parallel storage type
411 : #ifdef UG_PARALLEL
412 : A.set_storage_type(PST_ADDITIVE);
413 : A.set_layouts(dd->layouts());
414 : #endif
415 0 : }
416 :
417 : /**
418 : * This function adds the contributions of all passed element discretizations
419 : * on one given subset to the global Stiffness matrix.
420 : *
421 : * \param[in] vElemDisc element discretizations
422 : * \param[in] dd DoF Distribution
423 : * \param[in] si subset index
424 : * \param[in] bNonRegularGrid flag to indicate if non regular grid is used
425 : * \param[in,out] A Stiffness matrix
426 : * \param[in] u solution
427 : */
428 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
429 : template <typename TElem>
430 0 : void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::
431 : AssembleStiffnessMatrix( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
432 : ConstSmartPtr<DoFDistribution> dd,
433 : int si, bool bNonRegularGrid,
434 : matrix_type& A,
435 : const vector_type& u)
436 : {
437 : // check if only some elements are selected
438 0 : if(m_spAssTuner->selected_elements_used())
439 : {
440 : std::vector<TElem*> vElem;
441 0 : m_spAssTuner->collect_selected_elements(vElem, dd, si);
442 :
443 : // assembling is carried out only over those elements
444 : // which are selected and in subset si
445 : gass_type::template AssembleStiffnessMatrix<TElem>
446 0 : (vElemDisc, m_spApproxSpace->domain(), dd, vElem.begin(), vElem.end(), si,
447 : bNonRegularGrid, A, u, m_spAssTuner);
448 0 : }
449 : else
450 : {
451 : // general case: assembling over all elements in subset si
452 : gass_type::template AssembleStiffnessMatrix<TElem>
453 0 : (vElemDisc, m_spApproxSpace->domain(), dd,
454 : dd->template begin<TElem>(si), dd->template end<TElem>(si), si,
455 : bNonRegularGrid, A, u, m_spAssTuner);
456 : }
457 0 : }
458 :
459 : //////////////////////////////////////////////////////////////////////////////
460 : //////////////////////////////////////////////////////////////////////////////
461 : // Time Independent (stationary)
462 : //////////////////////////////////////////////////////////////////////////////
463 : //////////////////////////////////////////////////////////////////////////////
464 :
465 :
466 : ///////////////////////////////////////////////////////////////////////////////
467 : // Jacobian (stationary)
468 : ///////////////////////////////////////////////////////////////////////////////
469 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
470 0 : void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::
471 : assemble_jacobian(matrix_type& J,
472 : const vector_type& u,
473 : ConstSmartPtr<DoFDistribution> dd)
474 : {
475 : PROFILE_FUNC_GROUP("discretization");
476 : // update the elem discs
477 : update_disc_items();
478 0 : prep_assemble_loop(m_vElemDisc);
479 :
480 : // reset matrix to zero and resize
481 0 : m_spAssTuner->resize(dd, J);
482 :
483 : // Union of Subsets
484 0 : SubsetGroup unionSubsets;
485 : std::vector<SubsetGroup> vSSGrp;
486 :
487 : // pre process - modifies the solution, used for computing the defect
488 : const vector_type* pModifyU = &u;
489 : SmartPtr<vector_type> pModifyMemory;
490 0 : if( m_spAssTuner->modify_solution_enabled() ){
491 0 : pModifyMemory = u.clone();
492 : pModifyU = pModifyMemory.get();
493 : try{
494 0 : for(int type = 1; type < CT_ALL; type = type << 1){
495 0 : if(!(m_spAssTuner->constraint_type_enabled(type))) continue;
496 0 : for(size_t i = 0; i < m_vConstraint.size(); ++i)
497 0 : if(m_vConstraint[i]->type() & type)
498 0 : m_vConstraint[i]->modify_solution(*pModifyMemory, u, dd, type);
499 : }
500 0 : } UG_CATCH_THROW("Cannot modify solution.");
501 : }
502 :
503 :
504 :
505 : // create list of all subsets
506 : try{
507 0 : CreateSubsetGroups(vSSGrp, unionSubsets, m_vElemDisc, dd->subset_handler());
508 0 : }UG_CATCH_THROW("'DomainDiscretization': Can not create Subset Groups and Union.");
509 :
510 : // loop subsets
511 0 : for(size_t i = 0; i < unionSubsets.size(); ++i)
512 : {
513 : // get subset
514 : const int si = unionSubsets[i];
515 :
516 : // get dimension of the subset
517 0 : const int dim = DimensionOfSubset(*dd->subset_handler(), si);
518 :
519 : // request if subset is regular grid
520 0 : bool bNonRegularGrid = !unionSubsets.regular_grid(i);
521 :
522 : // overrule by regular grid if required
523 0 : if(m_spAssTuner->regular_grid_forced()) bNonRegularGrid = false;
524 :
525 : // Elem Disc on the subset
526 : std::vector<IElemDisc<TDomain>*> vSubsetElemDisc;
527 :
528 : // get all element discretizations that work on the subset
529 0 : GetElemDiscOnSubset(vSubsetElemDisc, m_vElemDisc, vSSGrp, si);
530 :
531 : // assemble on suitable elements
532 : try
533 : {
534 0 : switch(dim)
535 : {
536 0 : case 0:
537 : this->template AssembleJacobian<RegularVertex>
538 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, J, *pModifyU);
539 0 : break;
540 0 : case 1:
541 : this->template AssembleJacobian<RegularEdge>
542 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, J, *pModifyU);
543 : // When assembling over lower-dim manifolds that contain hanging nodes:
544 : this->template AssembleJacobian<ConstrainingEdge>
545 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, J, *pModifyU);
546 0 : break;
547 0 : case 2:
548 : this->template AssembleJacobian<Triangle>
549 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, J, *pModifyU);
550 : this->template AssembleJacobian<Quadrilateral>
551 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, J, *pModifyU);
552 : // When assembling over lower-dim manifolds that contain hanging nodes:
553 : this->template AssembleJacobian<ConstrainingTriangle>
554 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, J, *pModifyU);
555 : this->template AssembleJacobian<ConstrainingQuadrilateral>
556 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, J, *pModifyU);
557 0 : break;
558 0 : case 3:
559 : this->template AssembleJacobian<Tetrahedron>
560 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, J, *pModifyU);
561 : this->template AssembleJacobian<Pyramid>
562 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, J, *pModifyU);
563 : this->template AssembleJacobian<Prism>
564 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, J, *pModifyU);
565 : this->template AssembleJacobian<Hexahedron>
566 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, J, *pModifyU);
567 : this->template AssembleJacobian<Octahedron>
568 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, J, *pModifyU);
569 0 : break;
570 0 : default:
571 0 : UG_THROW("DomainDiscretization::assemble_jacobian (stationary):"
572 : "Dimension "<<dim<<"(subset="<<si<<") not supported");
573 : }
574 : }
575 0 : UG_CATCH_THROW("DomainDiscretization::assemble_jacobian (stationary):"
576 : " Assembling of elements of Dimension " << dim << " in "
577 : " subset "<<si<< " failed.");
578 : }
579 :
580 : // post process
581 : try{
582 0 : for(int type = 1; type < CT_ALL; type = type << 1){
583 0 : if(!(m_spAssTuner->constraint_type_enabled(type))) continue;
584 0 : for(size_t i = 0; i < m_vConstraint.size(); ++i)
585 0 : if(m_vConstraint[i]->type() & type)
586 : {
587 0 : m_vConstraint[i]->set_ass_tuner(m_spAssTuner);
588 0 : m_vConstraint[i]->adjust_jacobian(J, *pModifyU, dd, type);
589 : }
590 : }
591 0 : post_assemble_loop(m_vElemDisc);
592 0 : }UG_CATCH_THROW("DomainDiscretization::assemble_jacobian:"
593 : " Cannot execute post process.");
594 :
595 : // Remember parallel storage type
596 : #ifdef UG_PARALLEL
597 : J.set_storage_type(PST_ADDITIVE);
598 : J.set_layouts(dd->layouts());
599 : #endif
600 0 : }
601 :
602 : /**
603 : * This function adds the contributions of all passed element discretizations
604 : * on one given subset to the global Jacobian in the stationary case.
605 : *
606 : * \param[in] vElemDisc element discretizations
607 : * \param[in] si subset index
608 : * \param[in] bNonRegularGrid flag to indicate if non regular grid is used
609 : * \param[in,out] J jacobian
610 : * \param[in] u solution
611 : */
612 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
613 : template <typename TElem>
614 0 : void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::
615 : AssembleJacobian( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
616 : ConstSmartPtr<DoFDistribution> dd,
617 : int si, bool bNonRegularGrid,
618 : matrix_type& J,
619 : const vector_type& u)
620 : {
621 : // check if only some elements are selected
622 0 : if(m_spAssTuner->selected_elements_used())
623 : {
624 : std::vector<TElem*> vElem;
625 0 : m_spAssTuner->collect_selected_elements(vElem, dd, si);
626 :
627 : // assembling is carried out only over those elements
628 : // which are selected and in subset si
629 : gass_type::template AssembleJacobian<TElem>
630 0 : (vElemDisc, m_spApproxSpace->domain(), dd, vElem.begin(), vElem.end(), si,
631 : bNonRegularGrid, J, u, m_spAssTuner);
632 0 : }
633 : else
634 : {
635 : // general case: assembling over all elements in subset si
636 : gass_type::template AssembleJacobian<TElem>
637 0 : (vElemDisc, m_spApproxSpace->domain(), dd,
638 : dd->template begin<TElem>(si), dd->template end<TElem>(si), si,
639 : bNonRegularGrid, J, u, m_spAssTuner);
640 : }
641 0 : }
642 :
643 : ///////////////////////////////////////////////////////////////////////////////
644 : // Defect (stationary)
645 : ///////////////////////////////////////////////////////////////////////////////
646 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
647 0 : void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::
648 : assemble_defect(vector_type& d,
649 : const vector_type& u,
650 : ConstSmartPtr<DoFDistribution> dd)
651 : {
652 : PROFILE_FUNC_GROUP("discretization");
653 : // update the elem discs
654 : update_disc_items();
655 0 : prep_assemble_loop(m_vElemDisc);
656 :
657 : // reset matrix to zero and resize
658 0 : m_spAssTuner->resize(dd, d);
659 :
660 : // Union of Subsets
661 0 : SubsetGroup unionSubsets;
662 : std::vector<SubsetGroup> vSSGrp;
663 :
664 : // pre process - modifies the solution, used for computing the defect
665 : const vector_type* pModifyU = &u;
666 : SmartPtr<vector_type> pModifyMemory;
667 0 : if( m_spAssTuner->modify_solution_enabled() ){
668 0 : pModifyMemory = u.clone();
669 : pModifyU = pModifyMemory.get();
670 : try{
671 0 : for(int type = 1; type < CT_ALL; type = type << 1){
672 0 : if(!(m_spAssTuner->constraint_type_enabled(type))) continue;
673 0 : for(size_t i = 0; i < m_vConstraint.size(); ++i)
674 0 : if(m_vConstraint[i]->type() & type)
675 0 : m_vConstraint[i]->modify_solution(*pModifyMemory, u, dd, type);
676 : }
677 0 : } UG_CATCH_THROW("Cannot modify solution.");
678 : }
679 :
680 : // create list of all subsets
681 : try{
682 0 : CreateSubsetGroups(vSSGrp, unionSubsets, m_vElemDisc, dd->subset_handler());
683 0 : }UG_CATCH_THROW("'DomainDiscretization': Can not create Subset Groups and Union.");
684 :
685 : // loop subsets
686 0 : for(size_t i = 0; i < unionSubsets.size(); ++i)
687 : {
688 : // get subset
689 : const int si = unionSubsets[i];
690 :
691 : // get dimension of the subset
692 0 : const int dim = DimensionOfSubset(*dd->subset_handler(), si);
693 :
694 : // request if subset is regular grid
695 0 : bool bNonRegularGrid = !unionSubsets.regular_grid(i);
696 :
697 : // overrule by regular grid if required
698 0 : if(m_spAssTuner->regular_grid_forced()) bNonRegularGrid = false;
699 :
700 : // Elem Disc on the subset
701 : std::vector<IElemDisc<TDomain>*> vSubsetElemDisc;
702 :
703 : // get all element discretizations that work on the subset
704 0 : GetElemDiscOnSubset(vSubsetElemDisc, m_vElemDisc, vSSGrp, si);
705 :
706 : // assemble on suitable elements
707 : try
708 : {
709 0 : switch(dim)
710 : {
711 0 : case 0:
712 : this->template AssembleDefect<RegularVertex>
713 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, d, *pModifyU);
714 0 : break;
715 0 : case 1:
716 : this->template AssembleDefect<RegularEdge>
717 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, d, *pModifyU);
718 : // When assembling over lower-dim manifolds that contain hanging nodes:
719 : this->template AssembleDefect<ConstrainingEdge>
720 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, d, *pModifyU);
721 0 : break;
722 0 : case 2:
723 : this->template AssembleDefect<Triangle>
724 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, d, *pModifyU);
725 : this->template AssembleDefect<Quadrilateral>
726 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, d, *pModifyU);
727 : // When assembling over lower-dim manifolds that contain hanging nodes:
728 : this->template AssembleDefect<ConstrainingTriangle>
729 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, d, *pModifyU);
730 : this->template AssembleDefect<ConstrainingQuadrilateral>
731 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, d, *pModifyU);
732 0 : break;
733 0 : case 3:
734 : this->template AssembleDefect<Tetrahedron>
735 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, d, *pModifyU);
736 : this->template AssembleDefect<Pyramid>
737 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, d, *pModifyU);
738 : this->template AssembleDefect<Prism>
739 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, d, *pModifyU);
740 : this->template AssembleDefect<Hexahedron>
741 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, d, *pModifyU);
742 : this->template AssembleDefect<Octahedron>
743 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, d, *pModifyU);
744 0 : break;
745 0 : default:
746 0 : UG_THROW("DomainDiscretization::assemble_defect (stationary):"
747 : "Dimension "<<dim<<" (subset="<<si<<") not supported.");
748 : }
749 : }
750 0 : UG_CATCH_THROW("DomainDiscretization::assemble_defect (stationary):"
751 : " Assembling of elements of Dimension " << dim << " in "
752 : " subset "<<si<< " failed.");
753 : }
754 :
755 : // post process
756 : try{
757 :
758 : // Dirichlet first, since hanging nodes might be constrained by Dirichlet nodes
759 0 : if (m_spAssTuner->constraint_type_enabled(CT_DIRICHLET))
760 : {
761 0 : for (size_t i = 0; i < m_vConstraint.size(); ++i)
762 : {
763 0 : if (m_vConstraint[i]->type() & CT_DIRICHLET)
764 : {
765 0 : m_vConstraint[i]->set_ass_tuner(m_spAssTuner);
766 0 : m_vConstraint[i]->adjust_defect(d, *pModifyU, dd, CT_DIRICHLET);
767 : }
768 : }
769 : }
770 :
771 0 : for(int type = 1; type < CT_ALL; type = type << 1){
772 0 : if(!(m_spAssTuner->constraint_type_enabled(type))) continue;
773 0 : for(size_t i = 0; i < m_vConstraint.size(); ++i)
774 0 : if(m_vConstraint[i]->type() & type)
775 : {
776 0 : m_vConstraint[i]->set_ass_tuner(m_spAssTuner);
777 0 : m_vConstraint[i]->adjust_defect(d, *pModifyU, dd, type);
778 : }
779 : }
780 0 : post_assemble_loop(m_vElemDisc);
781 0 : } UG_CATCH_THROW("Cannot adjust defect.");
782 :
783 :
784 : // Remember parallel storage type
785 : #ifdef UG_PARALLEL
786 : d.set_storage_type(PST_ADDITIVE);
787 : #endif
788 0 : }
789 :
790 : /**
791 : * This function adds the contributions of all passed element discretizations
792 : * on one given subset to the global Defect in the stationary case.
793 : *
794 : * \param[in] vElemDisc element discretizations
795 : * \param[in] dd DoF Distribution
796 : * \param[in] si subset index
797 : * \param[in] bNonRegularGrid flag to indicate if non regular grid is used
798 : * \param[in,out] d defect
799 : * \param[in] u solution
800 : */
801 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
802 : template <typename TElem>
803 0 : void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::
804 : AssembleDefect( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
805 : ConstSmartPtr<DoFDistribution> dd,
806 : int si, bool bNonRegularGrid,
807 : vector_type& d,
808 : const vector_type& u)
809 : {
810 : // check if only some elements are selected
811 0 : if(m_spAssTuner->selected_elements_used())
812 : {
813 : std::vector<TElem*> vElem;
814 0 : m_spAssTuner->collect_selected_elements(vElem, dd, si);
815 :
816 : // assembling is carried out only over those elements
817 : // which are selected and in subset si
818 : gass_type::template AssembleDefect<TElem>
819 0 : (vElemDisc, m_spApproxSpace->domain(), dd, vElem.begin(), vElem.end(), si,
820 : bNonRegularGrid, d, u, m_spAssTuner);
821 0 : }
822 : else
823 : {
824 : // general case: assembling over all elements in subset si
825 : gass_type::template AssembleDefect<TElem>
826 0 : (vElemDisc, m_spApproxSpace->domain(), dd,
827 : dd->template begin<TElem>(si), dd->template end<TElem>(si), si,
828 : bNonRegularGrid, d, u, m_spAssTuner);
829 : }
830 0 : }
831 :
832 : ///////////////////////////////////////////////////////////////////////////////
833 : // Matrix and RHS (stationary)
834 : ///////////////////////////////////////////////////////////////////////////////
835 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
836 0 : void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::
837 : assemble_linear(matrix_type& mat, vector_type& rhs,
838 : ConstSmartPtr<DoFDistribution> dd)
839 : {
840 : PROFILE_FUNC_GROUP("discretization");
841 : // update the elem discs
842 : update_disc_items();
843 0 : prep_assemble_loop(m_vElemDisc);
844 :
845 : // reset matrix to zero and resize
846 0 : m_spAssTuner->resize(dd, mat);
847 0 : m_spAssTuner->resize(dd, rhs);
848 :
849 : // Union of Subsets
850 0 : SubsetGroup unionSubsets;
851 : std::vector<SubsetGroup> vSSGrp;
852 :
853 : // create list of all subsets
854 : try{
855 0 : CreateSubsetGroups(vSSGrp, unionSubsets, m_vElemDisc, dd->subset_handler());
856 0 : }UG_CATCH_THROW("DomainDiscretization: Can not create Subset Groups and Union.");
857 :
858 : // loop subsets
859 0 : for(size_t i = 0; i < unionSubsets.size(); ++i)
860 : {
861 : // get subset
862 : const int si = unionSubsets[i];
863 :
864 : // get dimension of the subset
865 0 : const int dim = DimensionOfSubset(*dd->subset_handler(), si);
866 :
867 : // request if subset is regular grid
868 0 : bool bNonRegularGrid = !unionSubsets.regular_grid(i);
869 :
870 : // overrule by regular grid if required
871 0 : if(m_spAssTuner->regular_grid_forced()) bNonRegularGrid = false;
872 :
873 : // Elem Disc on the subset
874 : std::vector<IElemDisc<TDomain>*> vSubsetElemDisc;
875 :
876 : // get all element discretizations that work on the subset
877 0 : GetElemDiscOnSubset(vSubsetElemDisc, m_vElemDisc, vSSGrp, si);
878 :
879 : // assemble on suitable elements
880 : try
881 : {
882 0 : switch(dim)
883 : {
884 0 : case 0:
885 : this->template AssembleLinear<RegularVertex>
886 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, mat, rhs);
887 0 : break;
888 0 : case 1:
889 : this->template AssembleLinear<RegularEdge>
890 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, mat, rhs);
891 : // When assembling over lower-dim manifolds that contain hanging nodes:
892 : this->template AssembleLinear<ConstrainingEdge>
893 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, mat, rhs);
894 0 : break;
895 0 : case 2:
896 : this->template AssembleLinear<Triangle>
897 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, mat, rhs);
898 : this->template AssembleLinear<Quadrilateral>
899 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, mat, rhs);
900 : // When assembling over lower-dim manifolds that contain hanging nodes:
901 : this->template AssembleLinear<ConstrainingTriangle>
902 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, mat, rhs);
903 : this->template AssembleLinear<ConstrainingQuadrilateral>
904 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, mat, rhs);
905 0 : break;
906 0 : case 3:
907 : this->template AssembleLinear<Tetrahedron>
908 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, mat, rhs);
909 : this->template AssembleLinear<Pyramid>
910 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, mat, rhs);
911 : this->template AssembleLinear<Prism>
912 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, mat, rhs);
913 : this->template AssembleLinear<Hexahedron>
914 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, mat, rhs);
915 : this->template AssembleLinear<Octahedron>
916 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, mat, rhs);
917 0 : break;
918 0 : default:
919 0 : UG_THROW("DomainDiscretization::assemble_linear (stationary):"
920 : "Dimension "<<dim<<" (subset="<<si<<") not supported.");
921 : }
922 : }
923 0 : UG_CATCH_THROW("DomainDiscretization::assemble_linear (stationary):"
924 : " Assembling of elements of Dimension " << dim << " in "
925 : " subset "<<si<< " failed.");
926 : }
927 :
928 : // post process
929 : try{
930 0 : for(int type = 1; type < CT_ALL; type = type << 1){
931 0 : if(!(m_spAssTuner->constraint_type_enabled(type))) continue;
932 0 : for(size_t i = 0; i < m_vConstraint.size(); ++i)
933 0 : if(m_vConstraint[i]->type() & type)
934 : {
935 0 : m_vConstraint[i]->set_ass_tuner(m_spAssTuner);
936 0 : m_vConstraint[i]->adjust_linear(mat, rhs, dd, type);
937 : }
938 : }
939 0 : post_assemble_loop(m_vElemDisc);
940 0 : }UG_CATCH_THROW("DomainDiscretization::assemble_linear: Cannot post process.");
941 :
942 : // Remember parallel storage type
943 : #ifdef UG_PARALLEL
944 : mat.set_storage_type(PST_ADDITIVE);
945 : mat.set_layouts(dd->layouts());
946 : rhs.set_storage_type(PST_ADDITIVE);
947 : #endif
948 0 : }
949 :
950 : /**
951 : * This function adds the contributions of all passed element discretizations
952 : * on one given subset to the global Matrix and the global Right-Hand Side
953 : * of the Linear problem in the stationary case.
954 : *
955 : * \param[in] vElemDisc element discretizations
956 : * \param[in] dd DoF Distribution
957 : * \param[in] si subset index
958 : * \param[in] bNonRegularGrid flag to indicate if non regular grid is used
959 : * \param[in,out] A Matrix
960 : * \param[in,out] rhs Right-hand side
961 : */
962 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
963 : template <typename TElem>
964 0 : void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::
965 : AssembleLinear( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
966 : ConstSmartPtr<DoFDistribution> dd,
967 : int si, bool bNonRegularGrid,
968 : matrix_type& A,
969 : vector_type& rhs)
970 : {
971 : // check if only some elements are selected
972 0 : if(m_spAssTuner->selected_elements_used())
973 : {
974 : std::vector<TElem*> vElem;
975 0 : m_spAssTuner->collect_selected_elements(vElem, dd, si);
976 :
977 : // assembling is carried out only over those elements
978 : // which are selected and in subset si
979 : gass_type::template AssembleLinear<TElem>
980 0 : (vElemDisc, m_spApproxSpace->domain(), dd, vElem.begin(), vElem.end(), si,
981 : bNonRegularGrid, A, rhs, m_spAssTuner);
982 0 : }
983 : else
984 : {
985 : // general case: assembling over all elements in subset si
986 : gass_type::template AssembleLinear<TElem>
987 0 : (vElemDisc, m_spApproxSpace->domain(), dd,
988 : dd->template begin<TElem>(si), dd->template end<TElem>(si), si,
989 : bNonRegularGrid, A, rhs, m_spAssTuner);
990 : }
991 0 : }
992 :
993 : ///////////////////////////////////////////////////////////////////////////////
994 : // RHS (stationary)
995 : ///////////////////////////////////////////////////////////////////////////////
996 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
997 0 : void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::
998 : assemble_rhs(vector_type& rhs,
999 : const vector_type& u,
1000 : ConstSmartPtr<DoFDistribution> dd)
1001 : {
1002 : PROFILE_FUNC_GROUP("discretization");
1003 : // update the elem discs
1004 : update_disc_items();
1005 0 : prep_assemble_loop(m_vElemDisc);
1006 :
1007 : // reset matrix to zero and resize
1008 0 : m_spAssTuner->resize(dd, rhs);
1009 :
1010 : // Union of Subsets
1011 0 : SubsetGroup unionSubsets;
1012 : std::vector<SubsetGroup> vSSGrp;
1013 :
1014 : // create list of all subsets
1015 : try{
1016 0 : CreateSubsetGroups(vSSGrp, unionSubsets, m_vElemDisc, dd->subset_handler());
1017 0 : }UG_CATCH_THROW("'DomainDiscretization': Can not create Subset Groups and Union.");
1018 :
1019 : // loop subsets
1020 0 : for(size_t i = 0; i < unionSubsets.size(); ++i)
1021 : {
1022 : // get subset
1023 : const int si = unionSubsets[i];
1024 :
1025 : // get dimension of the subset
1026 0 : const int dim = DimensionOfSubset(*dd->subset_handler(), si);
1027 :
1028 : // request if subset is regular grid
1029 0 : bool bNonRegularGrid = !unionSubsets.regular_grid(i);
1030 :
1031 : // overrule by regular grid if required
1032 0 : if(m_spAssTuner->regular_grid_forced()) bNonRegularGrid = false;
1033 :
1034 : // Elem Disc on the subset
1035 : std::vector<IElemDisc<TDomain>*> vSubsetElemDisc;
1036 :
1037 : // get all element discretizations that work on the subset
1038 0 : GetElemDiscOnSubset(vSubsetElemDisc, m_vElemDisc, vSSGrp, si);
1039 :
1040 : // assemble on suitable elements
1041 : try
1042 : {
1043 0 : switch(dim)
1044 : {
1045 0 : case 0:
1046 : this->template AssembleRhs<RegularVertex>
1047 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, rhs, u);
1048 0 : break;
1049 0 : case 1:
1050 : this->template AssembleRhs<RegularEdge>
1051 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, rhs, u);
1052 : // When assembling over lower-dim manifolds that contain hanging nodes:
1053 : this->template AssembleRhs<ConstrainingEdge>
1054 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, rhs, u);
1055 0 : break;
1056 0 : case 2:
1057 : this->template AssembleRhs<Triangle>
1058 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, rhs, u);
1059 : this->template AssembleRhs<Quadrilateral>
1060 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, rhs, u);
1061 : // When assembling over lower-dim manifolds that contain hanging nodes:
1062 : this->template AssembleRhs<ConstrainingTriangle>
1063 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, rhs, u);
1064 : this->template AssembleRhs<ConstrainingQuadrilateral>
1065 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, rhs, u);
1066 0 : break;
1067 0 : case 3:
1068 : this->template AssembleRhs<Tetrahedron>
1069 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, rhs, u);
1070 : this->template AssembleRhs<Pyramid>
1071 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, rhs, u);
1072 : this->template AssembleRhs<Prism>
1073 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, rhs, u);
1074 : this->template AssembleRhs<Hexahedron>
1075 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, rhs, u);
1076 : this->template AssembleRhs<Octahedron>
1077 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, rhs, u);
1078 0 : break;
1079 0 : default:
1080 0 : UG_THROW("DomainDiscretization::assemble_rhs (stationary):"
1081 : "Dimension "<<dim<<" (subset="<<si<<") not supported.");
1082 : }
1083 : }
1084 0 : UG_CATCH_THROW("DomainDiscretization::assemble_rhs (stationary):"
1085 : " Assembling of elements of Dimension " << dim << " in "
1086 : " subset "<<si<< " failed.");
1087 : }
1088 :
1089 : // post process
1090 : try{
1091 0 : for(int type = 1; type < CT_ALL; type = type << 1){
1092 0 : if(!(m_spAssTuner->constraint_type_enabled(type))) continue;
1093 0 : for(size_t i = 0; i < m_vConstraint.size(); ++i)
1094 0 : if(m_vConstraint[i]->type() & type)
1095 : {
1096 0 : m_vConstraint[i]->set_ass_tuner(m_spAssTuner);
1097 0 : m_vConstraint[i]->adjust_rhs(rhs, u, dd, type);
1098 : }
1099 : }
1100 0 : post_assemble_loop(m_vElemDisc);
1101 0 : }UG_CATCH_THROW("DomainDiscretization::assemble_rhs:"
1102 : " Cannot execute post process.");
1103 :
1104 : // Remember parallel storage type
1105 : #ifdef UG_PARALLEL
1106 : rhs.set_storage_type(PST_ADDITIVE);
1107 : #endif
1108 0 : }
1109 :
1110 : /**
1111 : * This function adds the contributions of all passed element discretizations
1112 : * on one given subset to the global Right-Hand Side.
1113 : *
1114 : * \param[in] vElemDisc element discretizations
1115 : * \param[in] dd DoF Distribution
1116 : * \param[in] si subset index
1117 : * \param[in] bNonRegularGrid flag to indicate if non regular grid is used
1118 : * \param[in,out] rhs Right-hand side
1119 : * \param[in] u solution
1120 : */
1121 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
1122 : template <typename TElem>
1123 0 : void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::
1124 : AssembleRhs( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
1125 : ConstSmartPtr<DoFDistribution> dd,
1126 : int si, bool bNonRegularGrid,
1127 : vector_type& rhs,
1128 : const vector_type& u)
1129 : {
1130 : // check if only some elements are selected
1131 0 : if(m_spAssTuner->selected_elements_used())
1132 : {
1133 : std::vector<TElem*> vElem;
1134 0 : m_spAssTuner->collect_selected_elements(vElem, dd, si);
1135 :
1136 : // assembling is carried out only over those elements
1137 : // which are selected and in subset si
1138 : gass_type::template AssembleRhs<TElem>
1139 0 : (vElemDisc, m_spApproxSpace->domain(), dd, vElem.begin(), vElem.end(), si,
1140 : bNonRegularGrid, rhs, u, m_spAssTuner);
1141 0 : }
1142 : else
1143 : {
1144 : // general case: assembling over all elements in subset si
1145 : gass_type::template AssembleRhs<TElem>
1146 0 : (vElemDisc, m_spApproxSpace->domain(), dd, dd->template begin<TElem>(si), dd->template end<TElem>(si), si,
1147 : bNonRegularGrid, rhs, u, m_spAssTuner);
1148 : }
1149 0 : }
1150 :
1151 : ///////////////////////////////////////////////////////////////////////////////
1152 : // set constraints (stationary)
1153 : ///////////////////////////////////////////////////////////////////////////////
1154 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
1155 0 : void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::
1156 : adjust_solution(vector_type& u, ConstSmartPtr<DoFDistribution> dd)
1157 : {
1158 : PROFILE_FUNC_GROUP("discretization");
1159 0 : update_constraints();
1160 :
1161 : // NOTE: it is crucial, that dirichlet pp are processed before constraints.
1162 : // otherwise we may start with an inconsistent solution in the solvers
1163 0 : std::vector<int> vType(6);
1164 0 : vType[0] = CT_DIRICHLET; // hanging (or other constrained) nodes might depend on Dirichlet nodes
1165 0 : vType[1] = CT_CONSTRAINTS;
1166 0 : vType[2] = CT_HANGING;
1167 0 : vType[3] = CT_MAY_DEPEND_ON_HANGING;
1168 0 : vType[4] = CT_ASSEMBLED;
1169 0 : vType[5] = CT_DIRICHLET; // hanging DoFs might be Dirichlet (Dirichlet overrides)
1170 :
1171 : // if assembling is carried out at one DoF only, u needs to be resized
1172 0 : if (m_spAssTuner->single_index_assembling_enabled()) u.resize(1);
1173 :
1174 : try{
1175 0 : for(size_t i = 0; i < vType.size(); ++i){
1176 0 : int type = vType[i];
1177 0 : if(!(m_spAssTuner->constraint_type_enabled(type))) continue;
1178 0 : for(size_t i = 0; i < m_vConstraint.size(); ++i)
1179 0 : if(m_vConstraint[i]->type() & type)
1180 : {
1181 0 : m_vConstraint[i]->set_ass_tuner(m_spAssTuner);
1182 0 : m_vConstraint[i]->adjust_solution(u, dd, type);
1183 : }
1184 : }
1185 :
1186 0 : } UG_CATCH_THROW("Cannot adjust solution.");
1187 0 : }
1188 :
1189 :
1190 :
1191 :
1192 : //////////////////////////////////////////////////////////////////////////////
1193 : //////////////////////////////////////////////////////////////////////////////
1194 : // Time Dependent (instationary)
1195 : //////////////////////////////////////////////////////////////////////////////
1196 : //////////////////////////////////////////////////////////////////////////////
1197 :
1198 : ///////////////////////////////////////////////////////////////////////////////
1199 : // Prepare Timestep (instationary)
1200 : ///////////////////////////////////////////////////////////////////////////////
1201 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
1202 0 : void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::
1203 : prepare_timestep
1204 : (
1205 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
1206 : number future_time,
1207 : ConstSmartPtr<DoFDistribution> dd
1208 : )
1209 : {
1210 : PROFILE_FUNC_GROUP("discretization");
1211 : // update the elem discs
1212 : update_disc_items();
1213 0 : prep_assemble_loop(m_vElemDisc);
1214 :
1215 : // find out whether grid is regular
1216 : ConstSmartPtr<ISubsetHandler> sh = dd->subset_handler();
1217 0 : size_t num_subsets = sh->num_subsets();
1218 : bool bNonRegularGrid = false;
1219 0 : for (size_t si = 0; si < num_subsets; ++si)
1220 0 : bNonRegularGrid |= !SubsetIsRegularGrid(*sh, si);
1221 :
1222 : // overrule by regular grid if required
1223 0 : if(m_spAssTuner->regular_grid_forced()) bNonRegularGrid = false;
1224 :
1225 : // call assembler's PrepareTimestep
1226 : try
1227 : {
1228 0 : gass_type::PrepareTimestep(m_vElemDisc, dd, bNonRegularGrid, vSol, future_time, m_spAssTuner);
1229 : }
1230 0 : UG_CATCH_THROW("DomainDiscretization::prepare_timestep (instationary):" <<
1231 : " Preparing time step failed.");
1232 :
1233 0 : post_assemble_loop(m_vElemDisc);
1234 0 : }
1235 :
1236 : ///////////////////////////////////////////////////////////////////////////////
1237 : // Prepare Timestep Elem (instationary)
1238 : ///////////////////////////////////////////////////////////////////////////////
1239 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
1240 0 : void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::
1241 : prepare_timestep_elem(ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
1242 : ConstSmartPtr<DoFDistribution> dd)
1243 : {
1244 : PROFILE_FUNC_GROUP("discretization");
1245 : // update the elem discs
1246 : update_disc_items();
1247 0 : prep_assemble_loop(m_vElemDisc);
1248 :
1249 : // Union of Subsets
1250 0 : SubsetGroup unionSubsets;
1251 : std::vector<SubsetGroup> vSSGrp;
1252 :
1253 : // create list of all subsets
1254 : try{
1255 0 : CreateSubsetGroups(vSSGrp, unionSubsets, m_vElemDisc, dd->subset_handler());
1256 0 : }UG_CATCH_THROW("'DomainDiscretization': Can not create Subset Groups and Union.");
1257 :
1258 : // loop subsets
1259 0 : for(size_t i = 0; i < unionSubsets.size(); ++i)
1260 : {
1261 : // get subset
1262 : const int si = unionSubsets[i];
1263 :
1264 : // get dimension of the subset
1265 0 : const int dim = DimensionOfSubset(*dd->subset_handler(), si);
1266 :
1267 : // request if subset is regular grid
1268 0 : bool bNonRegularGrid = !unionSubsets.regular_grid(i);
1269 :
1270 : // overrule by regular grid if required
1271 0 : if(m_spAssTuner->regular_grid_forced()) bNonRegularGrid = false;
1272 :
1273 : // Elem Disc on the subset
1274 : std::vector<IElemDisc<TDomain>*> vSubsetElemDisc;
1275 :
1276 : // get all element discretizations that work on the subset
1277 0 : GetElemDiscOnSubset(vSubsetElemDisc, m_vElemDisc, vSSGrp, si);
1278 :
1279 : // assemble on suitable elements
1280 : try
1281 : {
1282 0 : switch(dim)
1283 : {
1284 : case 0:
1285 : this->template PrepareTimestepElem<RegularVertex>
1286 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, vSol);
1287 0 : break;
1288 : case 1:
1289 : this->template PrepareTimestepElem<RegularEdge>
1290 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, vSol);
1291 : // When assembling over lower-dim manifolds that contain hanging nodes:
1292 : this->template PrepareTimestepElem<ConstrainingEdge>
1293 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, vSol);
1294 0 : break;
1295 : case 2:
1296 : this->template PrepareTimestepElem<Triangle>
1297 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, vSol);
1298 : this->template PrepareTimestepElem<Quadrilateral>
1299 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, vSol);
1300 : // When assembling over lower-dim manifolds that contain hanging nodes:
1301 : this->template PrepareTimestepElem<ConstrainingTriangle>
1302 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, vSol);
1303 : this->template PrepareTimestepElem<ConstrainingQuadrilateral>
1304 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, vSol);
1305 0 : break;
1306 : case 3:
1307 : this->template PrepareTimestepElem<Tetrahedron>
1308 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, vSol);
1309 : this->template PrepareTimestepElem<Pyramid>
1310 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, vSol);
1311 : this->template PrepareTimestepElem<Prism>
1312 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, vSol);
1313 : this->template PrepareTimestepElem<Hexahedron>
1314 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, vSol);
1315 : this->template PrepareTimestepElem<Octahedron>
1316 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, vSol);
1317 0 : break;
1318 0 : default:
1319 0 : UG_THROW("DomainDiscretization::prepare_timestep_elem (instationary):"
1320 : "Dimension "<<dim<<" (subset="<<si<<") not supported.");
1321 : }
1322 : }
1323 0 : UG_CATCH_THROW("DomainDiscretization::prepare_timestep_elem (instationary):"
1324 : " Assembling of elements of Dimension " << dim << " in "
1325 : " subset "<<si<< " failed.");
1326 : }
1327 0 : post_assemble_loop(m_vElemDisc);
1328 0 : }
1329 :
1330 : /**
1331 : * This function prepares the global discretization for a time-stepping scheme
1332 : * by calling the "prepare_timestep_elem" methods of all passed element
1333 : * discretizations on one given subset.
1334 : *
1335 : * \param[in] vElemDisc element discretizations
1336 : * \param[in] dd DoF Distribution
1337 : * \param[in] si subset index
1338 : * \param[in] bNonRegularGrid flag to indicate if non regular grid is used
1339 : * \param[in] vSol current and previous solutions
1340 : */
1341 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
1342 : template <typename TElem>
1343 0 : void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::
1344 : PrepareTimestepElem(const std::vector<IElemDisc<domain_type>*>& vElemDisc,
1345 : ConstSmartPtr<DoFDistribution> dd,
1346 : int si, bool bNonRegularGrid,
1347 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol)
1348 : {
1349 : // check if only some elements are selected
1350 0 : if(m_spAssTuner->selected_elements_used())
1351 : {
1352 : std::vector<TElem*> vElem;
1353 0 : m_spAssTuner->collect_selected_elements(vElem, dd, si);
1354 :
1355 : // assembling is carried out only over those elements
1356 : // which are selected and in subset si
1357 : gass_type::template PrepareTimestepElem<TElem>
1358 0 : (vElemDisc, m_spApproxSpace->domain(), dd, vElem.begin(), vElem.end(), si,
1359 : bNonRegularGrid, vSol, m_spAssTuner);
1360 0 : }
1361 : else
1362 : {
1363 : // general case: assembling over all elements in subset si
1364 : gass_type::template PrepareTimestepElem<TElem>
1365 0 : (vElemDisc, m_spApproxSpace->domain(), dd,
1366 : dd->template begin<TElem>(si), dd->template end<TElem>(si), si,
1367 : bNonRegularGrid, vSol, m_spAssTuner);
1368 : }
1369 0 : }
1370 :
1371 : ///////////////////////////////////////////////////////////////////////////////
1372 : // Jacobian (instationary)
1373 : ///////////////////////////////////////////////////////////////////////////////
1374 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
1375 0 : void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::
1376 : assemble_jacobian(matrix_type& J,
1377 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
1378 : const number s_a0,
1379 : ConstSmartPtr<DoFDistribution> dd)
1380 : {
1381 : // do not do anything if matrix is supposed to be const
1382 0 : if (m_spAssTuner->matrix_is_const()) return;
1383 :
1384 : PROFILE_FUNC_GROUP("discretization");
1385 : // update the elem discs
1386 : update_disc_items();
1387 0 : prep_assemble_loop(m_vElemDisc);
1388 :
1389 : // reset matrix to zero and resize
1390 0 : m_spAssTuner->resize(dd, J);
1391 :
1392 : // get current time
1393 : const number time = vSol->time(0);
1394 :
1395 : // Union of Subsets
1396 0 : SubsetGroup unionSubsets;
1397 : std::vector<SubsetGroup> vSSGrp;
1398 :
1399 : // create list of all subsets
1400 : try{
1401 0 : CreateSubsetGroups(vSSGrp, unionSubsets, m_vElemDisc, dd->subset_handler());
1402 0 : }UG_CATCH_THROW("'DomainDiscretization': Can not create Subset Groups and Union.");
1403 :
1404 : // preprocess - modifies the solution, used for computing the defect
1405 : ConstSmartPtr<VectorTimeSeries<vector_type> > pModifyU = vSol;
1406 : SmartPtr<VectorTimeSeries<vector_type> > pModifyMemory;
1407 0 : if( m_spAssTuner->modify_solution_enabled() ){
1408 0 : pModifyMemory = vSol->clone();
1409 0 : pModifyU = pModifyMemory;
1410 : try{
1411 0 : for(int type = 1; type < CT_ALL; type = type << 1){
1412 0 : if(!(m_spAssTuner->constraint_type_enabled(type))) continue;
1413 0 : for(size_t i = 0; i < m_vConstraint.size(); ++i)
1414 0 : if(m_vConstraint[i]->type() & type)
1415 0 : m_vConstraint[i]->modify_solution(pModifyMemory, vSol, dd, type);
1416 : }
1417 0 : } UG_CATCH_THROW("'DomainDiscretization': Cannot modify solution.");
1418 : }
1419 :
1420 : // loop subsets
1421 0 : for(size_t i = 0; i < unionSubsets.size(); ++i)
1422 : {
1423 : // get subset
1424 : const int si = unionSubsets[i];
1425 :
1426 : // get dimension of the subset
1427 0 : const int dim = DimensionOfSubset(*dd->subset_handler(), si);
1428 :
1429 : // request if subset is regular grid
1430 0 : bool bNonRegularGrid = !unionSubsets.regular_grid(i);
1431 :
1432 : // overrule by regular grid if required
1433 0 : if(m_spAssTuner->regular_grid_forced()) bNonRegularGrid = false;
1434 :
1435 : // Elem Disc on the subset
1436 : std::vector<IElemDisc<TDomain>*> vSubsetElemDisc;
1437 :
1438 : // get all element discretizations that work on the subset
1439 0 : GetElemDiscOnSubset(vSubsetElemDisc, m_vElemDisc, vSSGrp, si);
1440 :
1441 : // assemble on suitable elements
1442 : try
1443 : {
1444 0 : switch(dim)
1445 : {
1446 : case 0:
1447 : this->template AssembleJacobian<RegularVertex>
1448 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, J, pModifyU, s_a0);
1449 0 : break;
1450 : case 1:
1451 : this->template AssembleJacobian<RegularEdge>
1452 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, J, pModifyU, s_a0);
1453 : // When assembling over lower-dim manifolds that contain hanging nodes:
1454 : this->template AssembleJacobian<ConstrainingEdge>
1455 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, J, pModifyU, s_a0);
1456 0 : break;
1457 : case 2:
1458 : this->template AssembleJacobian<Triangle>
1459 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, J, pModifyU, s_a0);
1460 : this->template AssembleJacobian<Quadrilateral>
1461 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, J, pModifyU, s_a0);
1462 : // When assembling over lower-dim manifolds that contain hanging nodes:
1463 : this->template AssembleJacobian<ConstrainingTriangle>
1464 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, J, pModifyU, s_a0);
1465 : this->template AssembleJacobian<ConstrainingQuadrilateral>
1466 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, J, pModifyU, s_a0);
1467 0 : break;
1468 : case 3:
1469 : this->template AssembleJacobian<Tetrahedron>
1470 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, J, pModifyU, s_a0);
1471 : this->template AssembleJacobian<Pyramid>
1472 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, J, pModifyU, s_a0);
1473 : this->template AssembleJacobian<Prism>
1474 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, J, pModifyU, s_a0);
1475 : this->template AssembleJacobian<Hexahedron>
1476 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, J, pModifyU, s_a0);
1477 : this->template AssembleJacobian<Octahedron>
1478 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, J, pModifyU, s_a0);
1479 0 : break;
1480 0 : default:
1481 0 : UG_THROW("DomainDiscretization::assemble_jacobian (instationary):"
1482 : "Dimension "<<dim<<" (subset="<<si<<") not supported.");
1483 : }
1484 : }
1485 0 : UG_CATCH_THROW("DomainDiscretization::assemble_jacobian (instationary):"
1486 : " Assembling of elements of Dimension " << dim << " in "
1487 : " subset "<<si<< " failed.");
1488 : }
1489 :
1490 : // post process
1491 : try{
1492 0 : for(int type = 1; type < CT_ALL; type = type << 1){
1493 0 : if(!(m_spAssTuner->constraint_type_enabled(type))) continue;
1494 0 : for(size_t i = 0; i < m_vConstraint.size(); ++i)
1495 0 : if(m_vConstraint[i]->type() & type)
1496 : {
1497 0 : m_vConstraint[i]->set_ass_tuner(m_spAssTuner);
1498 0 : m_vConstraint[i]->adjust_jacobian(J, *pModifyU->solution(0), dd, type, time, pModifyU,s_a0);
1499 : }
1500 : }
1501 0 : post_assemble_loop(m_vElemDisc);
1502 0 : }UG_CATCH_THROW("Cannot adjust jacobian.");
1503 :
1504 : // Remember parallel storage type
1505 : #ifdef UG_PARALLEL
1506 : J.set_storage_type(PST_ADDITIVE);
1507 : J.set_layouts(dd->layouts());
1508 : #endif
1509 0 : }
1510 :
1511 : /**
1512 : * This function adds the contributions of all passed element discretizations
1513 : * on one given subset to the global Jacobian in the time-dependent case.
1514 : * Note, that s_m0 == 1
1515 : *
1516 : * \param[in] vElemDisc element discretizations
1517 : * \param[in] dd DoF Distribution
1518 : * \param[in] si subset index
1519 : * \param[in] bNonRegularGrid flag to indicate if non regular grid is used
1520 : * \param[in,out] J jacobian
1521 : * \param[in] vSol current and previous solutions
1522 : * \param[in] s_a0 scaling factor for stiffness part
1523 : */
1524 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
1525 : template <typename TElem>
1526 0 : void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::
1527 : AssembleJacobian( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
1528 : ConstSmartPtr<DoFDistribution> dd,
1529 : int si, bool bNonRegularGrid,
1530 : matrix_type& J,
1531 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
1532 : number s_a0)
1533 : {
1534 : // check if only some elements are selected
1535 0 : if(m_spAssTuner->selected_elements_used())
1536 : {
1537 : std::vector<TElem*> vElem;
1538 0 : m_spAssTuner->collect_selected_elements(vElem, dd, si);
1539 :
1540 : // assembling is carried out only over those elements
1541 : // which are selected and in subset si
1542 : gass_type::template AssembleJacobian<TElem>
1543 0 : (vElemDisc, m_spApproxSpace->domain(), dd, vElem.begin(), vElem.end(), si,
1544 : bNonRegularGrid, J, vSol, s_a0, m_spAssTuner);
1545 0 : }
1546 : else
1547 : {
1548 : gass_type::template AssembleJacobian<TElem>
1549 0 : (vElemDisc, m_spApproxSpace->domain(), dd,
1550 : dd->template begin<TElem>(si), dd->template end<TElem>(si), si,
1551 : bNonRegularGrid, J, vSol, s_a0, m_spAssTuner);
1552 : }
1553 0 : }
1554 :
1555 : ///////////////////////////////////////////////////////////////////////////////
1556 : // Defect (instationary)
1557 : ///////////////////////////////////////////////////////////////////////////////
1558 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
1559 0 : void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::
1560 : assemble_defect(vector_type& d,
1561 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
1562 : const std::vector<number>& vScaleMass,
1563 : const std::vector<number>& vScaleStiff,
1564 : ConstSmartPtr<DoFDistribution> dd)
1565 : {
1566 : PROFILE_FUNC_GROUP("discretization");
1567 : // update the elem discs
1568 : update_disc_items();
1569 0 : prep_assemble_loop(m_vElemDisc);
1570 :
1571 : // reset vector to zero and resize
1572 0 : m_spAssTuner->resize(dd, d);
1573 :
1574 : // Union of Subsets
1575 0 : SubsetGroup unionSubsets;
1576 : std::vector<SubsetGroup> vSSGrp;
1577 :
1578 : // create list of all subsets
1579 : try{
1580 0 : CreateSubsetGroups(vSSGrp, unionSubsets, m_vElemDisc, dd->subset_handler());
1581 0 : }UG_CATCH_THROW("'DomainDiscretization': Can not create Subset Groups and Union.");
1582 :
1583 : // pre process - modifies the solution, used for computing the defect
1584 : ConstSmartPtr<VectorTimeSeries<vector_type> > pModifyU = vSol;
1585 : SmartPtr<VectorTimeSeries<vector_type> > pModifyMemory;
1586 0 : if( m_spAssTuner->modify_solution_enabled() ){
1587 0 : pModifyMemory = vSol->clone();
1588 0 : pModifyU = pModifyMemory;
1589 : try{
1590 0 : for(int type = 1; type < CT_ALL; type = type << 1){
1591 0 : if(!(m_spAssTuner->constraint_type_enabled(type))) continue;
1592 0 : for(size_t i = 0; i < m_vConstraint.size(); ++i)
1593 0 : if(m_vConstraint[i]->type() & type)
1594 0 : m_vConstraint[i]->modify_solution(pModifyMemory, vSol, dd, type);
1595 : }
1596 0 : } UG_CATCH_THROW("'DomainDiscretization: Cannot modify solution.");
1597 : }
1598 :
1599 : // loop subsets
1600 0 : for(size_t i = 0; i < unionSubsets.size(); ++i)
1601 : {
1602 : // get subset
1603 : const int si = unionSubsets[i];
1604 :
1605 : // get dimension of the subset
1606 0 : const int dim = DimensionOfSubset(*dd->subset_handler(), si);
1607 :
1608 : // request if subset is regular grid
1609 0 : bool bNonRegularGrid = !unionSubsets.regular_grid(i);
1610 :
1611 : // overrule by regular grid if required
1612 0 : if(m_spAssTuner->regular_grid_forced()) bNonRegularGrid = false;
1613 :
1614 : // Elem Disc on the subset
1615 : std::vector<IElemDisc<TDomain>*> vSubsetElemDisc;
1616 :
1617 : // get all element discretizations that work on the subset
1618 0 : GetElemDiscOnSubset(vSubsetElemDisc, m_vElemDisc, vSSGrp, si);
1619 :
1620 : // assemble on suitable elements
1621 : try
1622 : {
1623 0 : switch(dim)
1624 : {
1625 : case 0:
1626 : this->template AssembleDefect<RegularVertex>
1627 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, d, pModifyU, vScaleMass, vScaleStiff);
1628 0 : break;
1629 : case 1:
1630 : this->template AssembleDefect<RegularEdge>
1631 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, d, pModifyU, vScaleMass, vScaleStiff);
1632 : // When assembling over lower-dim manifolds that contain hanging nodes:
1633 : this->template AssembleDefect<ConstrainingEdge>
1634 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, d, pModifyU, vScaleMass, vScaleStiff);
1635 0 : break;
1636 : case 2:
1637 : this->template AssembleDefect<Triangle>
1638 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, d, pModifyU, vScaleMass, vScaleStiff);
1639 : this->template AssembleDefect<Quadrilateral>
1640 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, d, pModifyU, vScaleMass, vScaleStiff);
1641 : // When assembling over lower-dim manifolds that contain hanging nodes:
1642 : this->template AssembleDefect<ConstrainingTriangle>
1643 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, d, pModifyU, vScaleMass, vScaleStiff);
1644 : this->template AssembleDefect<ConstrainingQuadrilateral>
1645 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, d, pModifyU, vScaleMass, vScaleStiff);
1646 0 : break;
1647 : case 3:
1648 : this->template AssembleDefect<Tetrahedron>
1649 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, d, pModifyU, vScaleMass, vScaleStiff);
1650 : this->template AssembleDefect<Pyramid>
1651 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, d, pModifyU, vScaleMass, vScaleStiff);
1652 : this->template AssembleDefect<Prism>
1653 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, d, pModifyU, vScaleMass, vScaleStiff);
1654 : this->template AssembleDefect<Hexahedron>
1655 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, d, pModifyU, vScaleMass, vScaleStiff);
1656 : this->template AssembleDefect<Octahedron>
1657 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, d, pModifyU, vScaleMass, vScaleStiff);
1658 0 : break;
1659 0 : default:
1660 0 : UG_THROW("DomainDiscretization::assemble_defect (instationary):"
1661 : "Dimension "<<dim<<" (subset="<<si<<") not supported.");
1662 : }
1663 : }
1664 0 : UG_CATCH_THROW("DomainDiscretization::assemble_defect (instationary):"
1665 : " Assembling of elements of Dimension " << dim << " in"
1666 : " subset "<< si << " failed.");
1667 : }
1668 :
1669 : // post process
1670 : try{
1671 0 : for(int type = 1; type < CT_ALL; type = type << 1){
1672 0 : if(!(m_spAssTuner->constraint_type_enabled(type))) continue;
1673 0 : for(size_t i = 0; i < m_vConstraint.size(); ++i)
1674 0 : if(m_vConstraint[i]->type() & type)
1675 : {
1676 0 : m_vConstraint[i]->set_ass_tuner(m_spAssTuner);
1677 0 : m_vConstraint[i]->adjust_defect(d, *pModifyU->solution(0), dd, type, pModifyU->time(0), pModifyU, &vScaleMass, &vScaleStiff);
1678 : }
1679 : }
1680 0 : post_assemble_loop(m_vElemDisc);
1681 0 : } UG_CATCH_THROW("Cannot adjust defect.");
1682 :
1683 : // Remember parallel storage type
1684 : #ifdef UG_PARALLEL
1685 : d.set_storage_type(PST_ADDITIVE);
1686 : #endif
1687 0 : }
1688 :
1689 : /*
1690 : * This function adds the contributions of all passed element discretizations
1691 : * on one given subset to the global Defect in the instationary case.
1692 : *
1693 : * \param[in] vElemDisc element discretizations
1694 : * \param[in] dd DoF Distribution
1695 : * \param[in] si subset index
1696 : * \param[in] bNonRegularGrid flag to indicate if non regular grid is used
1697 : * \param[in,out] d defect
1698 : * \param[in] vSol current and previous solutions
1699 : * \param[in] vScaleMass scaling factors for mass part
1700 : * \param[in] vScaleStiff scaling factors for stiffness part
1701 : */
1702 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
1703 : template <typename TElem>
1704 0 : void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::
1705 : AssembleDefect( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
1706 : ConstSmartPtr<DoFDistribution> dd,
1707 : int si, bool bNonRegularGrid,
1708 : vector_type& d,
1709 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
1710 : const std::vector<number>& vScaleMass,
1711 : const std::vector<number>& vScaleStiff)
1712 : {
1713 : // check if only some elements are selected
1714 0 : if(m_spAssTuner->selected_elements_used())
1715 : {
1716 : std::vector<TElem*> vElem;
1717 0 : m_spAssTuner->collect_selected_elements(vElem, dd, si);
1718 :
1719 : // assembling is carried out only over those elements
1720 : // which are selected and in subset si
1721 : gass_type::template AssembleDefect<TElem>
1722 0 : (vElemDisc, m_spApproxSpace->domain(), dd, vElem.begin(), vElem.end(), si,
1723 : bNonRegularGrid, d, vSol, vScaleMass, vScaleStiff, m_spAssTuner);
1724 0 : }
1725 : else
1726 : {
1727 : // general case: assembling over all elements in subset si
1728 : gass_type::template AssembleDefect<TElem>
1729 0 : (vElemDisc, m_spApproxSpace->domain(), dd,
1730 : dd->template begin<TElem>(si), dd->template end<TElem>(si), si,
1731 : bNonRegularGrid, d, vSol, vScaleMass, vScaleStiff, m_spAssTuner);
1732 : }
1733 0 : }
1734 :
1735 : ///////////////////////////////////////////////////////////////////////////////
1736 : // Matrix and RHS (instationary)
1737 : ///////////////////////////////////////////////////////////////////////////////
1738 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
1739 0 : void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::
1740 : assemble_linear(matrix_type& mat, vector_type& rhs,
1741 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
1742 : const std::vector<number>& vScaleMass,
1743 : const std::vector<number>& vScaleStiff,
1744 : ConstSmartPtr<DoFDistribution> dd)
1745 : {
1746 : PROFILE_FUNC_GROUP("discretization");
1747 : // update the elem discs
1748 : update_disc_items();
1749 :
1750 : // reset matrix to zero and resize
1751 0 : if (!m_spAssTuner->matrix_is_const())
1752 0 : m_spAssTuner->resize(dd, mat);
1753 0 : m_spAssTuner->resize(dd, rhs);
1754 :
1755 : // Union of Subsets
1756 0 : SubsetGroup unionSubsets;
1757 : std::vector<SubsetGroup> vSSGrp;
1758 :
1759 : // create list of all subsets
1760 : try{
1761 0 : CreateSubsetGroups(vSSGrp, unionSubsets, m_vElemDisc, dd->subset_handler());
1762 0 : }UG_CATCH_THROW("'DomainDiscretization': Can not create Subset Groups and Union.");
1763 :
1764 : // loop subsets
1765 0 : for(size_t i = 0; i < unionSubsets.size(); ++i)
1766 : {
1767 : // get subset
1768 : const int si = unionSubsets[i];
1769 :
1770 : // get dimension of the subset
1771 0 : const int dim = DimensionOfSubset(*dd->subset_handler(), si);
1772 :
1773 : // request if subset is regular grid
1774 0 : bool bNonRegularGrid = !unionSubsets.regular_grid(i);
1775 :
1776 : // overrule by regular grid if required
1777 0 : if(m_spAssTuner->regular_grid_forced()) bNonRegularGrid = false;
1778 :
1779 : // Elem Disc on the subset
1780 : std::vector<IElemDisc<TDomain>*> vSubsetElemDisc;
1781 :
1782 : // get all element discretizations that work on the subset
1783 0 : GetElemDiscOnSubset(vSubsetElemDisc, m_vElemDisc, vSSGrp, si);
1784 :
1785 : // assemble on suitable elements
1786 : try
1787 : {
1788 0 : switch(dim)
1789 : {
1790 : case 0:
1791 : this->template AssembleLinear<RegularVertex>
1792 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, mat, rhs, vSol, vScaleMass, vScaleStiff);
1793 0 : break;
1794 : case 1:
1795 : this->template AssembleLinear<RegularEdge>
1796 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, mat, rhs, vSol, vScaleMass, vScaleStiff);
1797 : // When assembling over lower-dim manifolds that contain hanging nodes:
1798 : this->template AssembleLinear<ConstrainingEdge>
1799 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, mat, rhs, vSol, vScaleMass, vScaleStiff);
1800 0 : break;
1801 : case 2:
1802 : this->template AssembleLinear<Triangle>
1803 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, mat, rhs, vSol, vScaleMass, vScaleStiff);
1804 : this->template AssembleLinear<Quadrilateral>
1805 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, mat, rhs, vSol, vScaleMass, vScaleStiff);
1806 : // When assembling over lower-dim manifolds that contain hanging nodes:
1807 : this->template AssembleLinear<ConstrainingTriangle>
1808 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, mat, rhs, vSol, vScaleMass, vScaleStiff);
1809 : this->template AssembleLinear<ConstrainingQuadrilateral>
1810 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, mat, rhs, vSol, vScaleMass, vScaleStiff);
1811 0 : break;
1812 : case 3:
1813 : this->template AssembleLinear<Tetrahedron>
1814 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, mat, rhs, vSol, vScaleMass, vScaleStiff);
1815 : this->template AssembleLinear<Pyramid>
1816 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, mat, rhs, vSol, vScaleMass, vScaleStiff);
1817 : this->template AssembleLinear<Prism>
1818 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, mat, rhs, vSol, vScaleMass, vScaleStiff);
1819 : this->template AssembleLinear<Hexahedron>
1820 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, mat, rhs, vSol, vScaleMass, vScaleStiff);
1821 : this->template AssembleLinear<Octahedron>
1822 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, mat, rhs, vSol, vScaleMass, vScaleStiff);
1823 0 : break;
1824 0 : default:
1825 0 : UG_THROW("DomainDiscretization::assemble_linear (instationary):"
1826 : "Dimension "<<dim<<" (subset="<<si<<") not supported.");
1827 : }
1828 : }
1829 0 : UG_CATCH_THROW("DomainDiscretization::assemble_linear (instationary):"
1830 : " Assembling of elements of Dimension " << dim << " in "
1831 : " subset "<<si<< " failed.");
1832 : }
1833 :
1834 : // post process
1835 : try{
1836 0 : for(int type = 1; type < CT_ALL; type = type << 1){
1837 0 : if(!(m_spAssTuner->constraint_type_enabled(type))) continue;
1838 0 : for(size_t i = 0; i < m_vConstraint.size(); ++i)
1839 0 : if(m_vConstraint[i]->type() & type)
1840 : {
1841 0 : m_vConstraint[i]->set_ass_tuner(m_spAssTuner);
1842 0 : m_vConstraint[i]->adjust_linear(mat, rhs, dd, type, vSol->time(0));
1843 : }
1844 : }
1845 0 : } UG_CATCH_THROW("Cannot adjust linear.");
1846 :
1847 : // Remember parallel storage type
1848 : #ifdef UG_PARALLEL
1849 : mat.set_storage_type(PST_ADDITIVE);
1850 : mat.set_layouts(dd->layouts());
1851 :
1852 : rhs.set_storage_type(PST_ADDITIVE);
1853 : #endif
1854 0 : }
1855 :
1856 : /**
1857 : * This function adds the contributions of all passed element discretizations
1858 : * on one given subset to the global Matrix and the global Right-Hand Side
1859 : * of the Linear problem in the stationary case.
1860 : *
1861 : * \param[in] vElemDisc element discretizations
1862 : * \param[in] dd DoF Distribution
1863 : * \param[in] si subset index
1864 : * \param[in] bNonRegularGrid flag to indicate if non regular grid is used
1865 : * \param[in,out] A Matrix
1866 : * \param[in,out] rhs Right-hand side
1867 : * \param[in] vSol current and previous solutions
1868 : * \param[in] vScaleMass scaling factors for mass part
1869 : * \param[in] vScaleStiff scaling factors for stiffness part
1870 : */
1871 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
1872 : template <typename TElem>
1873 0 : void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::
1874 : AssembleLinear( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
1875 : ConstSmartPtr<DoFDistribution> dd,
1876 : int si, bool bNonRegularGrid,
1877 : matrix_type& A,
1878 : vector_type& rhs,
1879 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
1880 : const std::vector<number>& vScaleMass,
1881 : const std::vector<number>& vScaleStiff)
1882 : {
1883 : // check if only some elements are selected
1884 0 : if(m_spAssTuner->selected_elements_used())
1885 : {
1886 : std::vector<TElem*> vElem;
1887 0 : m_spAssTuner->collect_selected_elements(vElem, dd, si);
1888 :
1889 : // assembling is carried out only over those elements
1890 : // which are selected and in subset si
1891 : gass_type::template AssembleLinear<TElem>
1892 0 : (vElemDisc, m_spApproxSpace->domain(), dd, vElem.begin(), vElem.end(), si,
1893 : bNonRegularGrid, A, rhs, vSol, vScaleMass, vScaleStiff, m_spAssTuner);
1894 0 : }
1895 : else
1896 : {
1897 : // general case: assembling over all elements in subset si
1898 : gass_type::template AssembleLinear<TElem>
1899 0 : (vElemDisc, m_spApproxSpace->domain(), dd,
1900 : dd->template begin<TElem>(si), dd->template end<TElem>(si), si,
1901 : bNonRegularGrid, A, rhs, vSol, vScaleMass, vScaleStiff, m_spAssTuner);
1902 : }
1903 0 : }
1904 :
1905 : ///////////////////////////////////////////////////////////////////////////////
1906 : // RHS (instationary)
1907 : ///////////////////////////////////////////////////////////////////////////////
1908 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
1909 0 : void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::
1910 : assemble_rhs(vector_type& rhs,
1911 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
1912 : const std::vector<number>& vScaleMass,
1913 : const std::vector<number>& vScaleStiff,
1914 : ConstSmartPtr<DoFDistribution> dd)
1915 : {
1916 : PROFILE_FUNC_GROUP("discretization");
1917 : // update the elem discs
1918 : update_disc_items();
1919 :
1920 : // reset vector to zero and resize
1921 0 : m_spAssTuner->resize(dd, rhs);
1922 :
1923 : // Union of Subsets
1924 0 : SubsetGroup unionSubsets;
1925 : std::vector<SubsetGroup> vSSGrp;
1926 :
1927 : // create list of all subsets
1928 : try{
1929 0 : CreateSubsetGroups(vSSGrp, unionSubsets, m_vElemDisc, dd->subset_handler());
1930 0 : }UG_CATCH_THROW("'DomainDiscretization': Can not create Subset Groups and Union.");
1931 :
1932 : // loop subsets
1933 0 : for(size_t i = 0; i < unionSubsets.size(); ++i)
1934 : {
1935 : // get subset
1936 : const int si = unionSubsets[i];
1937 :
1938 : // get dimension of the subset
1939 0 : const int dim = DimensionOfSubset(*dd->subset_handler(), si);
1940 :
1941 : // request if subset is regular grid
1942 0 : bool bNonRegularGrid = !unionSubsets.regular_grid(i);
1943 :
1944 : // overrule by regular grid if required
1945 0 : if(m_spAssTuner->regular_grid_forced()) bNonRegularGrid = false;
1946 :
1947 : // Elem Disc on the subset
1948 : std::vector<IElemDisc<TDomain>*> vSubsetElemDisc;
1949 :
1950 : // get all element discretizations that work on the subset
1951 0 : GetElemDiscOnSubset(vSubsetElemDisc, m_vElemDisc, vSSGrp, si);
1952 :
1953 : // assemble on suitable elements
1954 : try
1955 : {
1956 0 : switch(dim)
1957 : {
1958 : case 0:
1959 : this->template AssembleRhs<RegularVertex>
1960 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, rhs, vSol, vScaleMass, vScaleStiff);
1961 0 : break;
1962 : case 1:
1963 : this->template AssembleRhs<RegularEdge>
1964 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, rhs, vSol, vScaleMass, vScaleStiff);
1965 : // When assembling over lower-dim manifolds that contain hanging nodes:
1966 : this->template AssembleRhs<ConstrainingEdge>
1967 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, rhs, vSol, vScaleMass, vScaleStiff);
1968 0 : break;
1969 : case 2:
1970 : this->template AssembleRhs<Triangle>
1971 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, rhs, vSol, vScaleMass, vScaleStiff);
1972 : this->template AssembleRhs<Quadrilateral>
1973 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, rhs, vSol, vScaleMass, vScaleStiff);
1974 : // When assembling over lower-dim manifolds that contain hanging nodes:
1975 : this->template AssembleRhs<ConstrainingTriangle>
1976 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, rhs, vSol, vScaleMass, vScaleStiff);
1977 : this->template AssembleRhs<ConstrainingQuadrilateral>
1978 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, rhs, vSol, vScaleMass, vScaleStiff);
1979 0 : break;
1980 : case 3:
1981 : this->template AssembleRhs<Tetrahedron>
1982 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, rhs, vSol, vScaleMass, vScaleStiff);
1983 : this->template AssembleRhs<Pyramid>
1984 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, rhs, vSol, vScaleMass, vScaleStiff);
1985 : this->template AssembleRhs<Prism>
1986 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, rhs, vSol, vScaleMass, vScaleStiff);
1987 : this->template AssembleRhs<Hexahedron>
1988 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, rhs, vSol, vScaleMass, vScaleStiff);
1989 : this->template AssembleRhs<Octahedron>
1990 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, rhs, vSol, vScaleMass, vScaleStiff);
1991 0 : break;
1992 0 : default:
1993 0 : UG_THROW("DomainDiscretization::assemble_rhs (instationary):"
1994 : "Dimension "<<dim<<" (subset="<<si<<") not supported.");
1995 : }
1996 : }
1997 0 : UG_CATCH_THROW("DomainDiscretization::assemble_rhs (instationary):"
1998 : " Assembling of elements of Dimension " << dim << " in "
1999 : " subset "<<si<< " failed.");
2000 : }
2001 :
2002 :
2003 : // post process
2004 : try{
2005 0 : for(int type = 1; type < CT_ALL; type = type << 1){
2006 0 : if(!(m_spAssTuner->constraint_type_enabled(type))) continue;
2007 0 : for(size_t i = 0; i < m_vConstraint.size(); ++i)
2008 0 : if(m_vConstraint[i]->type() & type)
2009 : {
2010 0 : m_vConstraint[i]->set_ass_tuner(m_spAssTuner);
2011 0 : m_vConstraint[i]->adjust_rhs(rhs, *(vSol->solution(0)), dd, type, vSol->time(0));
2012 : }
2013 : }
2014 0 : } UG_CATCH_THROW("Cannot adjust linear.");
2015 :
2016 : // Remember parallel storage type
2017 : #ifdef UG_PARALLEL
2018 : rhs.set_storage_type(PST_ADDITIVE);
2019 : #endif
2020 0 : }
2021 :
2022 : /**
2023 : * This function adds the contributions of all passed element discretizations
2024 : * on one given subset to the global Right-Hand Side in the time-dependent case.
2025 : *
2026 : * \param[in] vElemDisc element discretizations
2027 : * \param[in] dd DoF Distribution
2028 : * \param[in] si subset index
2029 : * \param[in] bNonRegularGrid flag to indicate if non regular grid is used
2030 : * \param[in,out] rhs Right-hand side
2031 : * \param[in] vSol current and previous solutions
2032 : * \param[in] vScaleMass scaling factors for mass part
2033 : * \param[in] vScaleStiff scaling factors for stiffness part
2034 : */
2035 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
2036 : template <typename TElem>
2037 0 : void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::
2038 : AssembleRhs( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
2039 : ConstSmartPtr<DoFDistribution> dd,
2040 : int si, bool bNonRegularGrid,
2041 : vector_type& rhs,
2042 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
2043 : const std::vector<number>& vScaleMass,
2044 : const std::vector<number>& vScaleStiff)
2045 : {
2046 : // check if only some elements are selected
2047 0 : if(m_spAssTuner->selected_elements_used())
2048 : {
2049 : std::vector<TElem*> vElem;
2050 0 : m_spAssTuner->collect_selected_elements(vElem, dd, si);
2051 :
2052 : // assembling is carried out only over those elements
2053 : // which are selected and in subset si
2054 : gass_type::template AssembleRhs<TElem>
2055 0 : (vElemDisc, m_spApproxSpace->domain(), dd, vElem.begin(), vElem.end(), si,
2056 : bNonRegularGrid, rhs, vSol, vScaleMass, vScaleStiff, m_spAssTuner);
2057 0 : }
2058 : else
2059 : {
2060 : // general case: assembling over all elements in subset si
2061 : gass_type::template AssembleRhs<TElem>
2062 0 : (vElemDisc, m_spApproxSpace->domain(), dd,
2063 : dd->template begin<TElem>(si), dd->template end<TElem>(si), si,
2064 : bNonRegularGrid, rhs, vSol, vScaleMass, vScaleStiff, m_spAssTuner);
2065 : }
2066 0 : }
2067 :
2068 : ///////////////////////////////////////////////////////////////////////////////
2069 : // set constraint values (instationary)
2070 : ///////////////////////////////////////////////////////////////////////////////
2071 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
2072 0 : void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::
2073 : adjust_solution(vector_type& u, number time, ConstSmartPtr<DoFDistribution> dd)
2074 : {
2075 : PROFILE_FUNC_GROUP("discretization");
2076 0 : update_constraints();
2077 :
2078 : // NOTE: it is crucial, that dirichlet pp are processed before constraints.
2079 : // otherwise we may start with an inconsistent solution in the solvers
2080 0 : std::vector<int> vType(5);
2081 0 : vType[0] = CT_DIRICHLET; // hanging (or other constrained) nodes might depend on Dirichlet nodes
2082 0 : vType[1] = CT_CONSTRAINTS;
2083 0 : vType[2] = CT_HANGING;
2084 0 : vType[3] = CT_MAY_DEPEND_ON_HANGING;
2085 0 : vType[4] = CT_DIRICHLET; // hanging DoFs might be Dirichlet (Dirichlet overrides)
2086 :
2087 : // if assembling is carried out at one DoF only, u needs to be resized
2088 0 : if (m_spAssTuner->single_index_assembling_enabled()) u.resize(1);
2089 :
2090 : try{
2091 :
2092 : // constraints
2093 0 : for(size_t i = 0; i < vType.size(); ++i){
2094 0 : int type = vType[i];
2095 0 : if(!(m_spAssTuner->constraint_type_enabled(type))) continue;
2096 0 : for(size_t i = 0; i < m_vConstraint.size(); ++i)
2097 0 : if(m_vConstraint[i]->type() & type)
2098 : {
2099 0 : m_vConstraint[i]->set_ass_tuner(m_spAssTuner);
2100 0 : m_vConstraint[i]->adjust_solution(u, dd, type, time);
2101 : }
2102 : }
2103 0 : } UG_CATCH_THROW(" Cannot adjust solution.");
2104 0 : }
2105 :
2106 : ///////////////////////////////////////////////////////////////////////////////
2107 : // Initialization of the exports (optional)
2108 : ///////////////////////////////////////////////////////////////////////////////
2109 :
2110 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
2111 0 : void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::
2112 : init_all_exports(ConstSmartPtr<DoFDistribution> dd, bool bAsTimeDependent)
2113 : {
2114 : PROFILE_FUNC_GROUP("discretization");
2115 : // update the elem discs
2116 : update_disc_items();
2117 :
2118 : // Union of Subsets
2119 0 : SubsetGroup unionSubsets;
2120 : std::vector<SubsetGroup> vSSGrp;
2121 :
2122 : // create list of all subsets
2123 : try{
2124 0 : CreateSubsetGroups(vSSGrp, unionSubsets, m_vElemDisc, dd->subset_handler());
2125 0 : }UG_CATCH_THROW("'DomainDiscretization': Can not create Subset Groups and Union.");
2126 :
2127 : // loop subsets
2128 0 : for(size_t i = 0; i < unionSubsets.size(); ++i)
2129 : {
2130 : // get subset
2131 : const int si = unionSubsets[i];
2132 :
2133 : // get dimension of the subset
2134 0 : const int dim = DimensionOfSubset(*dd->subset_handler(), si);
2135 :
2136 : // request if subset is regular grid
2137 0 : bool bNonRegularGrid = !unionSubsets.regular_grid(i);
2138 :
2139 : // overrule by regular grid if required
2140 0 : if(m_spAssTuner->regular_grid_forced()) bNonRegularGrid = false;
2141 :
2142 : // Elem Disc on the subset
2143 : std::vector<IElemDisc<TDomain>*> vSubsetElemDisc;
2144 :
2145 : // get all element discretizations that work on the subset
2146 0 : GetElemDiscOnSubset(vSubsetElemDisc, m_vElemDisc, vSSGrp, si);
2147 :
2148 : // assemble on suitable elements
2149 : try
2150 : {
2151 0 : switch(dim)
2152 : {
2153 0 : case 0:
2154 : this->template InitAllExports<RegularVertex>
2155 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, bAsTimeDependent);
2156 0 : break;
2157 0 : case 1:
2158 : this->template InitAllExports<RegularEdge>
2159 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, bAsTimeDependent);
2160 : // When assembling over lower-dim manifolds that contain hanging nodes:
2161 : this->template InitAllExports<ConstrainingEdge>
2162 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, bAsTimeDependent);
2163 0 : break;
2164 0 : case 2:
2165 : this->template InitAllExports<Triangle>
2166 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, bAsTimeDependent);
2167 : this->template InitAllExports<Quadrilateral>
2168 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, bAsTimeDependent);
2169 : // When assembling over lower-dim manifolds that contain hanging nodes:
2170 : this->template InitAllExports<ConstrainingTriangle>
2171 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, bAsTimeDependent);
2172 : this->template InitAllExports<ConstrainingQuadrilateral>
2173 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, bAsTimeDependent);
2174 0 : break;
2175 0 : case 3:
2176 : this->template InitAllExports<Tetrahedron>
2177 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, bAsTimeDependent);
2178 : this->template InitAllExports<Pyramid>
2179 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, bAsTimeDependent);
2180 : this->template InitAllExports<Prism>
2181 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, bAsTimeDependent);
2182 : this->template InitAllExports<Hexahedron>
2183 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, bAsTimeDependent);
2184 : this->template InitAllExports<Octahedron>
2185 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, bAsTimeDependent);
2186 0 : break;
2187 0 : default:
2188 0 : UG_THROW("DomainDiscretization::init_all_exports:"
2189 : "Dimension "<<dim<<" (subset="<<si<<") not supported.");
2190 : }
2191 : }
2192 0 : UG_CATCH_THROW("DomainDiscretization::assemble_stiffness_matrix:"
2193 : " Assembling of elements of Dimension " << dim << " in "
2194 : " subset "<<si<< " failed.");
2195 : }
2196 :
2197 0 : }
2198 :
2199 : /**
2200 : * This function initalizes the exports and does not do anything else.
2201 : *
2202 : * It transfers the dof distributions to the exports and creates the function
2203 : * maps there.
2204 : *
2205 : * This initizalization is optional but after that, the exports can be used
2206 : * by other objects (like in the plotting).
2207 : *
2208 : * \param[in] vElemDisc element discretizations
2209 : * \param[in] dd DoF Distribution
2210 : * \param[in] si subset index
2211 : * \param[in] bNonRegularGrid flag to indicate if non regular grid is used
2212 : * \param[in] bAsTimeDependent flag to simulate the instationary case for the initialization
2213 : */
2214 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
2215 : template <typename TElem>
2216 0 : void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::
2217 : InitAllExports( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
2218 : ConstSmartPtr<DoFDistribution> dd,
2219 : int si, bool bNonRegularGrid, bool bAsTimeDependent)
2220 : {
2221 : // check if only some elements are selected
2222 0 : if(m_spAssTuner->selected_elements_used())
2223 : {
2224 : std::vector<TElem*> vElem;
2225 0 : m_spAssTuner->collect_selected_elements(vElem, dd, si);
2226 :
2227 : // assembling is carried out only over those elements
2228 : // which are selected and in subset si
2229 : gass_type::template InitAllExports<TElem>
2230 0 : (vElemDisc, dd, vElem.begin(), vElem.end(), si, bNonRegularGrid, bAsTimeDependent);
2231 0 : }
2232 : else
2233 : {
2234 : // general case: assembling over all elements in subset si
2235 : gass_type::template InitAllExports<TElem>
2236 0 : (vElemDisc, dd,
2237 : dd->template begin<TElem>(si), dd->template end<TElem>(si), si,
2238 : bNonRegularGrid, bAsTimeDependent);
2239 : }
2240 0 : }
2241 :
2242 : ///////////////////////////////////////////////////////////////////////////////
2243 : ///////////////////////////////////////////////////////////////////////////////
2244 : // Error estimator
2245 : ///////////////////////////////////////////////////////////////////////////////
2246 : ///////////////////////////////////////////////////////////////////////////////
2247 :
2248 : template <typename TDomain, typename T>
2249 0 : void CollectIErrEstData(std::vector<IErrEstData<TDomain>*> &vErrEstData, const T &vElemDisc)
2250 : {
2251 0 : for (std::size_t i = 0; i < vElemDisc.size(); ++i)
2252 : {
2253 0 : SmartPtr<IErrEstData<TDomain> > sp_err_est_data = vElemDisc[i]->err_est_data();
2254 0 : IErrEstData<TDomain>* err_est_data = sp_err_est_data.get();
2255 0 : if (err_est_data == NULL) continue; // no data specified
2256 0 : if (std::find (vErrEstData.begin(), vErrEstData.end(), err_est_data) != vErrEstData.end())
2257 0 : continue; // this one is already in the array
2258 0 : if (err_est_data->consider_me()) vErrEstData.push_back(err_est_data);
2259 : }
2260 :
2261 0 : }
2262 :
2263 : ///////////////////////////////////////////////////////////////////////////////
2264 : // Error estimator (stationary)
2265 : ///////////////////////////////////////////////////////////////////////////////
2266 :
2267 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
2268 0 : void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::
2269 : calc_error
2270 : (
2271 : const vector_type& u,
2272 : ConstSmartPtr<DoFDistribution> dd,
2273 : error_vector_type* u_vtk
2274 : )
2275 : {
2276 : PROFILE_FUNC_GROUP("error_estimator");
2277 :
2278 : // get multigrid
2279 : SmartPtr<MultiGrid> pMG = ((DoFDistribution *) dd.get())->multi_grid();
2280 :
2281 : // check, whether separate error data exists
2282 : const bool useErrorData = !m_vDomainElemError.empty();
2283 : // UG_LOG("useErrorData=" << (useErrorData) << std::endl);
2284 :
2285 : // update the elem discs
2286 0 : if (useErrorData) { update_error_items();}
2287 : else { update_disc_items();}
2288 :
2289 : // Union of Subsets
2290 0 : SubsetGroup unionSubsets;
2291 : std::vector<SubsetGroup> vSSGrp;
2292 :
2293 : // create list of all subsets
2294 : try
2295 : {
2296 0 : if (useErrorData) { CreateSubsetGroups(vSSGrp, unionSubsets, m_vElemError, dd->subset_handler());}
2297 0 : else { CreateSubsetGroups(vSSGrp, unionSubsets, m_vElemDisc, dd->subset_handler());}
2298 : }
2299 0 : UG_CATCH_THROW("'DomainDiscretization': Can not create Subset Groups and Union.");
2300 :
2301 : // get the error estimator data for all the discretizations
2302 : std::vector<IErrEstData<TDomain>*> vErrEstData;
2303 0 : if (useErrorData) CollectIErrEstData(vErrEstData, m_vElemError);
2304 0 : else CollectIErrEstData(vErrEstData, m_vElemDisc);
2305 : /*for(std::size_t i = 0; i < m_vElemDisc.size(); ++i)
2306 : {
2307 : SmartPtr<IErrEstData<TDomain> > sp_err_est_data = m_vElemDisc[i]->err_est_data();
2308 : IErrEstData<TDomain>* err_est_data = sp_err_est_data.get();
2309 : if (err_est_data == NULL) continue; // no data specified
2310 : if (std::find (vErrEstData.begin(), vErrEstData.end(), err_est_data) != vErrEstData.end())
2311 : continue; // this one is already in the array
2312 : if (err_est_data->consider_me()) vErrEstData.push_back(err_est_data);
2313 : }*/
2314 :
2315 : // preprocess the error estimator data in the discretizations
2316 : try
2317 : {
2318 0 : for (size_t i = 0; i < vErrEstData.size(); ++i)
2319 0 : vErrEstData[i]->alloc_err_est_data(dd->surface_view(), dd->grid_level());
2320 : }
2321 0 : UG_CATCH_THROW("DomainDiscretization::calc_error: Cannot prepare the error estimator");
2322 :
2323 :
2324 :
2325 : // loop subsets to assemble the estimators
2326 0 : for (size_t i = 0; i < unionSubsets.size(); ++i)
2327 : {
2328 : // get subset
2329 : const int si = unionSubsets[i];
2330 :
2331 : // get dimension of the subset
2332 0 : const int dim = DimensionOfSubset(*dd->subset_handler(), si);
2333 :
2334 : // request if subset is regular grid
2335 0 : bool bNonRegularGrid = !unionSubsets.regular_grid(i);
2336 :
2337 : // Elem Disc on the subset
2338 : typedef typename std::vector<IElemError<TDomain>*> error_vector_type;
2339 : error_vector_type vSubsetElemError;
2340 :
2341 : // get all element discretizations that work on the subset
2342 0 : GetElemDiscItemOnSubset<IElemError<TDomain>, IElemDisc<TDomain> > (vSubsetElemError, m_vElemDisc, vSSGrp, si);
2343 0 : if (useErrorData)
2344 : {
2345 : GetElemDiscItemOnSubset<IElemError<TDomain>, IElemError<TDomain> >
2346 0 : (vSubsetElemError, m_vElemError, vSSGrp, si);
2347 : }
2348 : else
2349 : {
2350 : GetElemDiscItemOnSubset<IElemError<TDomain>, IElemDisc<TDomain> >
2351 0 : (vSubsetElemError, m_vElemDisc, vSSGrp, si);
2352 : }
2353 :
2354 : /*
2355 : UG_LOG("m_vElemDisc.size="<<m_vElemDisc.size()<< std::endl);
2356 : UG_LOG("m_vElemError.size="<<m_vElemError.size()<< std::endl);
2357 : UG_LOG("vSubsetElemError.size="<<vSubsetElemError.size()<< std::endl);
2358 : */
2359 :
2360 : // remove from elemDisc list those with !err_est_enabled()
2361 : typename error_vector_type::iterator it = vSubsetElemError.begin();
2362 0 : while (it != vSubsetElemError.end())
2363 : {
2364 0 : if (!(*it)->err_est_enabled())
2365 : it = vSubsetElemError.erase(it);
2366 : else ++it;
2367 : }
2368 : // UG_LOG("vSubsetElemError.size="<<vSubsetElemError.size()<< std::endl);
2369 :
2370 : // assemble on suitable elements
2371 : try
2372 : {
2373 0 : switch (dim)
2374 : {
2375 0 : case 1:
2376 : this->template AssembleErrorEstimator<RegularEdge>
2377 0 : (vSubsetElemError, dd, si, bNonRegularGrid, u);
2378 : // When assembling over lower-dim manifolds that contain hanging nodes:
2379 : this->template AssembleErrorEstimator<ConstrainingEdge>
2380 0 : (vSubsetElemError, dd, si, bNonRegularGrid, u);
2381 0 : break;
2382 0 : case 2:
2383 : this->template AssembleErrorEstimator<Triangle>
2384 0 : (vSubsetElemError, dd, si, bNonRegularGrid, u);
2385 : this->template AssembleErrorEstimator<Quadrilateral>
2386 0 : (vSubsetElemError, dd, si, bNonRegularGrid, u);
2387 : // When assembling over lower-dim manifolds that contain hanging nodes:
2388 : this->template AssembleErrorEstimator<ConstrainingTriangle>
2389 0 : (vSubsetElemError, dd, si, bNonRegularGrid, u);
2390 : this->template AssembleErrorEstimator<ConstrainingQuadrilateral>
2391 0 : (vSubsetElemError, dd, si, bNonRegularGrid, u);
2392 0 : break;
2393 0 : case 3:
2394 : this->template AssembleErrorEstimator<Tetrahedron>
2395 0 : (vSubsetElemError, dd, si, bNonRegularGrid, u);
2396 : this->template AssembleErrorEstimator<Pyramid>
2397 0 : (vSubsetElemError, dd, si, bNonRegularGrid, u);
2398 : this->template AssembleErrorEstimator<Prism>
2399 0 : (vSubsetElemError, dd, si, bNonRegularGrid, u);
2400 : this->template AssembleErrorEstimator<Hexahedron>
2401 0 : (vSubsetElemError, dd, si, bNonRegularGrid, u);
2402 : this->template AssembleErrorEstimator<Octahedron>
2403 0 : (vSubsetElemError, dd, si, bNonRegularGrid, u);
2404 0 : break;
2405 0 : default:
2406 0 : UG_THROW("DomainDiscretization::calc_error:"
2407 : " Dimension "<<dim<<" (subset="<<si<<") not supported.");
2408 : }
2409 : }
2410 0 : UG_CATCH_THROW("DomainDiscretization::calc_error:"
2411 : " Assembling of elements of Dimension " << dim << " in "
2412 : " subset "<<si<< " failed.");
2413 : }
2414 :
2415 : // apply constraints
2416 : try
2417 : {
2418 0 : for (int type = 1; type < CT_ALL; type = type << 1)
2419 : {
2420 0 : if (!(m_spAssTuner->constraint_type_enabled(type))) continue;
2421 0 : for (size_t i = 0; i < m_vConstraint.size(); ++i)
2422 0 : if ((m_vConstraint[i]->type() & type) && m_vConstraint[i]->err_est_enabled())
2423 : {
2424 0 : m_vConstraint[i]->set_ass_tuner(m_spAssTuner);
2425 0 : m_vConstraint[i]->adjust_error(u, dd, type);
2426 : }
2427 : }
2428 : }
2429 0 : UG_CATCH_THROW("Cannot adjust error assemblings.");
2430 :
2431 : // summarize the error estimator data in the discretizations
2432 : try
2433 : {
2434 0 : for (std::size_t i = 0; i < vErrEstData.size(); ++i)
2435 0 : vErrEstData[i]->summarize_err_est_data(m_spApproxSpace->domain());
2436 : }
2437 0 : UG_CATCH_THROW("DomainDiscretization::calc_error: Cannot summarize the error estimator");
2438 :
2439 : // perform integrations for error estimators and mark elements
2440 : typedef typename domain_traits<dim>::element_type elem_type;
2441 : typedef typename SurfaceView::traits<elem_type>::const_iterator elem_iter_type;
2442 :
2443 0 : m_mgElemErrors.attach_indicators(pMG);
2444 :
2445 : // loop surface elements
2446 : ConstSmartPtr<SurfaceView> sv = dd->surface_view();
2447 : const GridLevel& gl = dd->grid_level();
2448 : elem_iter_type elem_iter_end = sv->template end<elem_type> (gl, SurfaceView::ALL);
2449 0 : for (elem_iter_type elem = sv->template begin<elem_type> (gl, SurfaceView::ALL); elem != elem_iter_end; ++elem)
2450 : {
2451 : // clear attachment (to be on the safe side)
2452 0 : m_mgElemErrors.error(*elem) = 0.0;
2453 :
2454 : // get corner coordinates
2455 0 : std::vector<MathVector<dim> > vCornerCoords = std::vector<MathVector<dim> >(0);
2456 0 : CollectCornerCoordinates(vCornerCoords, *elem, m_spApproxSpace->domain()->position_accessor(), false);
2457 :
2458 : // integrate for all estimators, then add up
2459 0 : for (std::size_t ee = 0; ee < vErrEstData.size(); ++ee)
2460 0 : m_mgElemErrors.error(*elem) += vErrEstData[ee]->scaling_factor()
2461 0 : * vErrEstData[ee]->get_elem_error_indicator(*elem, &vCornerCoords[0]);
2462 : }
2463 :
2464 : // write error estimator values to vtk
2465 0 : if (u_vtk)
2466 : {
2467 : // local indices and solution
2468 : LocalIndices ind; LocalVector locU;
2469 :
2470 : // cast u_vtk to grid_function
2471 0 : GridFunction<TDomain, CPUAlgebra>* uVTK = dynamic_cast<GridFunction<TDomain, CPUAlgebra>*>(u_vtk);
2472 0 : if (!uVTK)
2473 : {
2474 0 : UG_THROW("Argument passed as output for error function is not a GridFunction of suitable type.");
2475 : }
2476 :
2477 : // clear previous values
2478 : uVTK->set(0.0);
2479 :
2480 : // map attachments to grid function
2481 0 : ConstSmartPtr<SurfaceView> sv = uVTK->approx_space()->dof_distribution(gl)->surface_view();
2482 : elem_iter_type elem_iter_end = sv->template end<elem_type> (gl, SurfaceView::ALL);
2483 0 : for (elem_iter_type elem = sv->template begin<elem_type> (gl, SurfaceView::ALL); elem != elem_iter_end; ++elem)
2484 : {
2485 : // get global indices
2486 0 : uVTK->approx_space()->dof_distribution(gl)->indices(*elem, ind, false);
2487 :
2488 0 : UG_COND_THROW(ind.num_fct() != 1,
2489 : "Number of functions in grid function passed for error indicator values is not 1 on "
2490 : << ElementDebugInfo(*uVTK->domain()->grid(), *elem) << ".");
2491 :
2492 0 : UG_COND_THROW(ind.num_dof(0) != 1,
2493 : "Number of DoFs in grid function passed for error indicator values is not 1 on "
2494 : << ElementDebugInfo(*uVTK->domain()->grid(), *elem) << ".");
2495 :
2496 : // adapt local algebra
2497 0 : locU.resize(ind);
2498 :
2499 : // read local values of u
2500 0 : GetLocalVector(locU, *uVTK);
2501 :
2502 : // assign error value
2503 0 : locU(0,0) = m_mgElemErrors.error(*elem);
2504 :
2505 : // add to grid function
2506 0 : AddLocalVector(*uVTK, locU);
2507 : }
2508 : }
2509 :
2510 0 : this->m_bErrorCalculated = true;
2511 :
2512 : // postprocess the error estimators in the discretizations
2513 : try{
2514 0 : for(std::size_t i = 0; i < vErrEstData.size(); ++i)
2515 0 : vErrEstData[i]->release_err_est_data();
2516 : }
2517 0 : UG_CATCH_THROW("DomainDiscretization::calc_error: Cannot release the error estimator");
2518 0 : }
2519 :
2520 : /**
2521 : * This function assembles the error estimator associated with all the
2522 : * element discretizations in the internal data structure for one given subset.
2523 : *
2524 : * \param[in] vElemDisc element discretizations
2525 : * \param[in] dd DoF Distribution
2526 : * \param[in] si subset index
2527 : * \param[in] bNonRegularGrid flag to indicate if non regular grid is used
2528 : * \param[in] u solution
2529 : */
2530 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
2531 : template <typename TElem>
2532 0 : inline void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::
2533 : AssembleErrorEstimator( const std::vector<IElemError<domain_type>*>& vElemDisc,
2534 : ConstSmartPtr<DoFDistribution> dd,
2535 : int si, bool bNonRegularGrid,
2536 : const vector_type& u)
2537 : {
2538 : // general case: assembling over all elements in subset si
2539 : gass_type::template AssembleErrorEstimator<TElem>
2540 0 : (vElemDisc, m_spApproxSpace->domain(), dd,
2541 : dd->template begin<TElem>(si), dd->template end<TElem>(si),
2542 : si, bNonRegularGrid, u);
2543 0 : }
2544 :
2545 : ///////////////////////////////////////////////////////////////////////////////
2546 : // Error estimator (instationary)
2547 : ///////////////////////////////////////////////////////////////////////////////
2548 :
2549 :
2550 :
2551 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
2552 0 : void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::
2553 : calc_error(ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
2554 : ConstSmartPtr<DoFDistribution> dd,
2555 : const std::vector<number>& vScaleMass,
2556 : const std::vector<number>& vScaleStiff,
2557 : error_vector_type* u_vtk)
2558 : {
2559 : PROFILE_FUNC_GROUP("error_estimator");
2560 :
2561 : // get multigrid
2562 : SmartPtr<MultiGrid> pMG = ((DoFDistribution *) dd.get())->multi_grid();
2563 :
2564 : // check, whether separate error data exists
2565 : const bool useErrorData = !m_vDomainElemError.empty();
2566 :
2567 : // update elem items
2568 0 : if (useErrorData) { update_error_items();}
2569 : else { update_disc_items();}
2570 :
2571 : // Union of Subsets
2572 0 : SubsetGroup unionSubsets;
2573 : std::vector<SubsetGroup> vSSGrp;
2574 :
2575 : // create list of all subsets
2576 : try
2577 : {
2578 0 : if (useErrorData) {CreateSubsetGroups(vSSGrp, unionSubsets, m_vElemError, dd->subset_handler());}
2579 0 : else {CreateSubsetGroups(vSSGrp, unionSubsets, m_vElemDisc, dd->subset_handler());}
2580 : }
2581 0 : UG_CATCH_THROW("'DomainDiscretization': Can not create Subset Groups and Union.");
2582 :
2583 : // collect the error estimator data (for all discretizations)
2584 : std::vector<IErrEstData<TDomain>*> vErrEstData;
2585 :
2586 0 : if (useErrorData) CollectIErrEstData(vErrEstData, m_vElemError);
2587 0 : else CollectIErrEstData(vErrEstData, m_vElemDisc);
2588 :
2589 : // preprocess the error estimator data in the discretizations
2590 : try
2591 : {
2592 0 : for (std::size_t i = 0; i < vErrEstData.size(); ++i)
2593 0 : vErrEstData[i]->alloc_err_est_data(dd->surface_view(), dd->grid_level());
2594 : }
2595 0 : UG_CATCH_THROW("DomainDiscretization::calc_error: Cannot prepare the error estimator");
2596 :
2597 : // loop subsets to assemble the estimators
2598 0 : for (std::size_t i = 0; i < unionSubsets.size(); ++i)
2599 : {
2600 : // get subset
2601 : const int si = unionSubsets[i];
2602 :
2603 : // get dimension of the subset
2604 0 : const int dim = DimensionOfSubset(*dd->subset_handler(), si);
2605 :
2606 : // request if subset is regular grid
2607 0 : bool bNonRegularGrid = !unionSubsets.regular_grid(i);
2608 :
2609 : // collect all element discretizations that work on the subset
2610 : typedef std::vector<IElemError<TDomain>*> error_vector_type;
2611 : error_vector_type vSubsetElemError;
2612 :
2613 0 : if (useErrorData) {
2614 : GetElemDiscItemOnSubset<IElemError<TDomain>, IElemError<TDomain> >
2615 0 : (vSubsetElemError, m_vElemError, vSSGrp, si);
2616 : }
2617 : else
2618 : {
2619 : GetElemDiscItemOnSubset<IElemError<TDomain>, IElemDisc<TDomain> >
2620 0 : (vSubsetElemError, m_vElemDisc, vSSGrp, si);
2621 : }
2622 :
2623 : // remove from elemDisc list those with !err_est_enabled()
2624 : typename error_vector_type::iterator it = vSubsetElemError.begin();
2625 0 : while (it != vSubsetElemError.end())
2626 : {
2627 0 : if (!(*it)->err_est_enabled())
2628 : it = vSubsetElemError.erase(it);
2629 : else ++it;
2630 : }
2631 :
2632 : // assemble on suitable elements
2633 : try
2634 : {
2635 0 : switch (dim)
2636 : {
2637 : case 1:
2638 : this->template AssembleErrorEstimator<RegularEdge>
2639 0 : (vSubsetElemError, dd, si, bNonRegularGrid, vScaleMass, vScaleStiff, vSol);
2640 : // When assembling over lower-dim manifolds that contain hanging nodes:
2641 : this->template AssembleErrorEstimator<ConstrainingEdge>
2642 0 : (vSubsetElemError, dd, si, bNonRegularGrid, vScaleMass, vScaleStiff, vSol);
2643 0 : break;
2644 : case 2:
2645 : this->template AssembleErrorEstimator<Triangle>
2646 0 : (vSubsetElemError, dd, si, bNonRegularGrid, vScaleMass, vScaleStiff, vSol);
2647 : this->template AssembleErrorEstimator<Quadrilateral>
2648 0 : (vSubsetElemError, dd, si, bNonRegularGrid, vScaleMass, vScaleStiff, vSol);
2649 : // When assembling over lower-dim manifolds that contain hanging nodes:
2650 : this->template AssembleErrorEstimator<ConstrainingTriangle>
2651 0 : (vSubsetElemError, dd, si, bNonRegularGrid, vScaleMass, vScaleStiff, vSol);
2652 : this->template AssembleErrorEstimator<ConstrainingQuadrilateral>
2653 0 : (vSubsetElemError, dd, si, bNonRegularGrid, vScaleMass, vScaleStiff, vSol);
2654 0 : break;
2655 : case 3:
2656 : this->template AssembleErrorEstimator<Tetrahedron>
2657 0 : (vSubsetElemError, dd, si, bNonRegularGrid, vScaleMass, vScaleStiff, vSol);
2658 : this->template AssembleErrorEstimator<Pyramid>
2659 0 : (vSubsetElemError, dd, si, bNonRegularGrid, vScaleMass, vScaleStiff, vSol);
2660 : this->template AssembleErrorEstimator<Prism>
2661 0 : (vSubsetElemError, dd, si, bNonRegularGrid, vScaleMass, vScaleStiff, vSol);
2662 : this->template AssembleErrorEstimator<Hexahedron>
2663 0 : (vSubsetElemError, dd, si, bNonRegularGrid, vScaleMass, vScaleStiff, vSol);
2664 : this->template AssembleErrorEstimator<Octahedron>
2665 0 : (vSubsetElemError, dd, si, bNonRegularGrid, vScaleMass, vScaleStiff, vSol);
2666 0 : break;
2667 0 : default:
2668 0 : UG_THROW("DomainDiscretization::calc_error:"
2669 : " Dimension "<<dim<<" (subset="<<si<<") not supported.");
2670 : }
2671 : }
2672 0 : UG_CATCH_THROW("DomainDiscretization::calc_error:"
2673 : " Assembling of elements of Dimension " << dim << " in "
2674 : " subset "<< si << "failed.");
2675 : }
2676 :
2677 : // apply constraints
2678 : try
2679 : {
2680 0 : for (int type = 1; type < CT_ALL; type = type << 1)
2681 : {
2682 0 : if (!(m_spAssTuner->constraint_type_enabled(type))) continue;
2683 0 : for (size_t i = 0; i < m_vConstraint.size(); ++i)
2684 0 : if ((m_vConstraint[i]->type() & type) && m_vConstraint[i]->err_est_enabled())
2685 : {
2686 0 : m_vConstraint[i]->set_ass_tuner(m_spAssTuner);
2687 0 : m_vConstraint[i]->adjust_error(*vSol->solution(0), dd, type, vSol->time(0),
2688 : vSol, &vScaleMass, &vScaleStiff);
2689 : }
2690 : }
2691 : }
2692 0 : UG_CATCH_THROW("Cannot adjust error assemblings.");
2693 :
2694 : // summarize the error estimator data in the discretizations
2695 : try
2696 : {
2697 0 : for (std::size_t i = 0; i < vErrEstData.size(); ++i)
2698 0 : vErrEstData[i]->summarize_err_est_data(m_spApproxSpace->domain());
2699 : }
2700 0 : UG_CATCH_THROW("DomainDiscretization::calc_error: Cannot summarize the error estimator.");
2701 :
2702 : // perform integrations for error estimators and mark elements
2703 : typedef typename domain_traits<dim>::element_type elem_type;
2704 : typedef typename SurfaceView::traits<elem_type>::const_iterator elem_iter_type;
2705 :
2706 0 : m_mgElemErrors.attach_indicators(pMG);
2707 :
2708 : // loop surface elements
2709 : ConstSmartPtr<SurfaceView> sv = dd->surface_view();
2710 : const GridLevel& gl = dd->grid_level();
2711 : elem_iter_type elem_iter_end = sv->template end<elem_type> (gl, SurfaceView::ALL);
2712 0 : for (elem_iter_type elem = sv->template begin<elem_type> (gl, SurfaceView::ALL); elem != elem_iter_end; ++elem)
2713 : {
2714 : // clear attachment
2715 0 : m_mgElemErrors.error(*elem) = 0.0;
2716 :
2717 : // get corner coordinates
2718 0 : std::vector<MathVector<dim> > vCornerCoords = std::vector<MathVector<dim> >(0);
2719 0 : CollectCornerCoordinates(vCornerCoords, *elem, m_spApproxSpace->domain()->position_accessor(), false);
2720 :
2721 : // integrate for all estimators, then add up
2722 0 : for (std::size_t ee = 0; ee < vErrEstData.size(); ++ee)
2723 0 : m_mgElemErrors.error(*elem) += vErrEstData[ee]->scaling_factor()
2724 0 : * vErrEstData[ee]->get_elem_error_indicator(*elem, &vCornerCoords[0]);
2725 : }
2726 :
2727 : // write error estimator values to vtk
2728 0 : if (u_vtk)
2729 : {
2730 : // local indices and solution
2731 : LocalIndices ind; LocalVector locU;
2732 :
2733 : // cast u_vtk to grid_function
2734 0 : GridFunction<TDomain,TAlgebra>* uVTK = dynamic_cast<GridFunction<TDomain,TAlgebra>*>(u_vtk);
2735 0 : if (!uVTK)
2736 : {
2737 0 : UG_THROW("Argument passed as output for error function is not a GridFunction.");
2738 : }
2739 :
2740 : // clear previous values
2741 : uVTK->set(0.0);
2742 :
2743 : // map attachments to grid function
2744 0 : ConstSmartPtr<SurfaceView> sv = uVTK->approx_space()->dof_distribution(gl)->surface_view();
2745 : elem_iter_type elem_iter_end = sv->template end<elem_type> (gl, SurfaceView::ALL);
2746 0 : for (elem_iter_type elem = sv->template begin<elem_type> (gl, SurfaceView::ALL); elem != elem_iter_end; ++elem)
2747 : {
2748 : // get global indices
2749 0 : uVTK->approx_space()->dof_distribution(gl)->indices(*elem, ind, false);
2750 :
2751 : // adapt local algebra
2752 0 : locU.resize(ind);
2753 :
2754 : // read local values of u
2755 0 : GetLocalVector(locU, *uVTK);
2756 :
2757 : // assign error value
2758 0 : locU(0,0) = m_mgElemErrors.error(*elem);
2759 :
2760 : // add to grid function
2761 0 : AddLocalVector(*uVTK, locU);
2762 : }
2763 : }
2764 :
2765 0 : this->m_bErrorCalculated = true;
2766 :
2767 : // postprocess the error estimators in the discretizations
2768 : try{
2769 0 : for(std::size_t i = 0; i < vErrEstData.size(); ++i)
2770 0 : vErrEstData[i]->release_err_est_data();
2771 : }
2772 0 : UG_CATCH_THROW("DomainDiscretization::calc_error: Cannot release the error estimator");
2773 0 : }
2774 :
2775 : /**
2776 : * This function assembles the error estimator associated with all the
2777 : * element discretizations in the internal data structure for one given subset.
2778 : *
2779 : * \param[in] vElemDisc element discretizations
2780 : * \param[in] dd DoF Distribution
2781 : * \param[in] si subset index
2782 : * \param[in] bNonRegularGrid flag to indicate if non regular grid is used
2783 : * \param[in] vScaleMass scaling factors for mass part
2784 : * \param[in] vScaleStiff scaling factors for stiffness part
2785 : * \param[in] vSol solution
2786 : */
2787 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
2788 : template <typename TElem>
2789 0 : inline void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::
2790 : AssembleErrorEstimator( const std::vector<IElemError<domain_type>*>& vElemDisc,
2791 : ConstSmartPtr<DoFDistribution> dd,
2792 : int si, bool bNonRegularGrid,
2793 : const std::vector<number>& vScaleMass,
2794 : const std::vector<number>& vScaleStiff,
2795 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol)
2796 : {
2797 : // general case: assembling over all elements in subset si
2798 : gass_type::template AssembleErrorEstimator<TElem>
2799 0 : (vElemDisc, m_spApproxSpace->domain(), dd,
2800 : dd->template begin<TElem>(si), dd->template end<TElem>(si),
2801 : si, bNonRegularGrid, vScaleMass, vScaleStiff, vSol);
2802 0 : }
2803 :
2804 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
2805 0 : void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::
2806 : mark_with_strategy
2807 : ( IRefiner& refiner,
2808 : SmartPtr<IElementMarkingStrategy<TDomain> >spMarkingStrategy
2809 : )
2810 : {
2811 : // check that error indicators have been calculated
2812 0 : if (!this->m_bErrorCalculated)
2813 : {
2814 0 : UG_THROW("Error indicators have to be calculated first by a call to 'calc_error'.");
2815 : }
2816 :
2817 : // mark elements for refinement
2818 0 : if (spMarkingStrategy.valid())
2819 : {
2820 0 : spMarkingStrategy->mark( m_mgElemErrors, refiner, this->dd(GridLevel(GridLevel::TOP, GridLevel::SURFACE)));
2821 : }
2822 0 : }
2823 :
2824 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
2825 0 : void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::
2826 : invalidate_error()
2827 : {
2828 : typedef typename domain_traits<dim>::element_type elem_type;
2829 :
2830 : // check that error indicators have been calculated
2831 0 : if (m_bErrorCalculated)
2832 : {
2833 0 : m_bErrorCalculated = false;
2834 0 : m_mgElemErrors.detach_indicators();
2835 : // this->m_pMG->template detach_from<elem_type>(this->m_aError);
2836 : }
2837 0 : }
2838 :
2839 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
2840 0 : bool DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::
2841 : is_error_valid()
2842 : {
2843 : // check that error indicators have been calculated
2844 0 : return this->m_bErrorCalculated;
2845 : }
2846 :
2847 : ///////////////////////////////////////////////////////////////////////////////
2848 : // Finish Timestep (instationary)
2849 : ///////////////////////////////////////////////////////////////////////////////
2850 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
2851 0 : void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::
2852 : finish_timestep
2853 : (
2854 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
2855 : ConstSmartPtr<DoFDistribution> dd
2856 : )
2857 : {
2858 : PROFILE_FUNC_GROUP("discretization");
2859 : // update the elem discs
2860 : update_disc_items();
2861 0 : prep_assemble_loop(m_vElemDisc);
2862 :
2863 : // find out whether grid is regular
2864 : ConstSmartPtr<ISubsetHandler> sh = dd->subset_handler();
2865 0 : size_t num_subsets = sh->num_subsets();
2866 : bool bNonRegularGrid = false;
2867 0 : for (size_t si = 0; si < num_subsets; ++si)
2868 0 : bNonRegularGrid |= !SubsetIsRegularGrid(*sh, si);
2869 :
2870 : // overrule by regular grid if required
2871 0 : if(m_spAssTuner->regular_grid_forced()) bNonRegularGrid = false;
2872 :
2873 : // call assembler's FinishTimestep
2874 : try
2875 : {
2876 0 : gass_type::FinishTimestep(m_vElemDisc, dd, bNonRegularGrid, vSol, m_spAssTuner);
2877 : }
2878 0 : UG_CATCH_THROW("DomainDiscretization::finish_timestep (instationary):" <<
2879 : " Finishing time step failed.");
2880 :
2881 0 : post_assemble_loop(m_vElemDisc);
2882 0 : }
2883 :
2884 : ///////////////////////////////////////////////////////////////////////////////
2885 : // Finish Timestep Elem (instationary)
2886 : ///////////////////////////////////////////////////////////////////////////////
2887 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
2888 0 : void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::
2889 : finish_timestep_elem(ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
2890 : ConstSmartPtr<DoFDistribution> dd)
2891 : {
2892 : PROFILE_FUNC_GROUP("discretization");
2893 : // update the elem discs
2894 : update_disc_items();
2895 :
2896 : // Union of Subsets
2897 0 : SubsetGroup unionSubsets;
2898 : std::vector<SubsetGroup> vSSGrp;
2899 :
2900 : // create list of all subsets
2901 : try{
2902 0 : CreateSubsetGroups(vSSGrp, unionSubsets, m_vElemDisc, dd->subset_handler());
2903 0 : }UG_CATCH_THROW("'DomainDiscretization': Can not create Subset Groups and Union.");
2904 :
2905 : // loop subsets
2906 0 : for(size_t i = 0; i < unionSubsets.size(); ++i)
2907 : {
2908 : // get subset
2909 : const int si = unionSubsets[i];
2910 :
2911 : // get dimension of the subset
2912 0 : const int dim = DimensionOfSubset(*dd->subset_handler(), si);
2913 :
2914 : // request if subset is regular grid
2915 0 : bool bNonRegularGrid = !unionSubsets.regular_grid(i);
2916 :
2917 : // overrule by regular grid if required
2918 0 : if(m_spAssTuner->regular_grid_forced()) bNonRegularGrid = false;
2919 :
2920 : // Elem Disc on the subset
2921 : std::vector<IElemDisc<TDomain>*> vSubsetElemDisc;
2922 :
2923 : // get all element discretizations that work on the subset
2924 0 : GetElemDiscOnSubset(vSubsetElemDisc, m_vElemDisc, vSSGrp, si);
2925 :
2926 : // assemble on suitable elements
2927 : try
2928 : {
2929 0 : switch(dim)
2930 : {
2931 : case 0:
2932 : this->template FinishTimestepElem<RegularVertex>
2933 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, vSol);
2934 0 : break;
2935 : case 1:
2936 : this->template FinishTimestepElem<RegularEdge>
2937 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, vSol);
2938 : // When assembling over lower-dim manifolds that contain hanging nodes:
2939 : this->template FinishTimestepElem<ConstrainingEdge>
2940 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, vSol);
2941 0 : break;
2942 : case 2:
2943 : this->template FinishTimestepElem<Triangle>
2944 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, vSol);
2945 : this->template FinishTimestepElem<Quadrilateral>
2946 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, vSol);
2947 : // When assembling over lower-dim manifolds that contain hanging nodes:
2948 : this->template FinishTimestepElem<ConstrainingTriangle>
2949 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, vSol);
2950 : this->template FinishTimestepElem<ConstrainingQuadrilateral>
2951 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, vSol);
2952 0 : break;
2953 : case 3:
2954 : this->template FinishTimestepElem<Tetrahedron>
2955 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, vSol);
2956 : this->template FinishTimestepElem<Pyramid>
2957 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, vSol);
2958 : this->template FinishTimestepElem<Prism>
2959 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, vSol);
2960 : this->template FinishTimestepElem<Hexahedron>
2961 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, vSol);
2962 : this->template FinishTimestepElem<Octahedron>
2963 0 : (vSubsetElemDisc, dd, si, bNonRegularGrid, vSol);
2964 0 : break;
2965 0 : default:
2966 0 : UG_THROW("DomainDiscretization::finish_timestep_elem (instationary):"
2967 : "Dimension "<<dim<<" (subset="<<si<<") not supported.");
2968 : }
2969 : }
2970 0 : UG_CATCH_THROW("DomainDiscretization::finish_timestep_elem (instationary):"
2971 : " Assembling of elements of Dimension " << dim << " in "
2972 : " subset "<<si<< " failed.");
2973 : }
2974 :
2975 0 : }
2976 :
2977 : /**
2978 : * This function finalizes the global discretization in a time-stepping scheme
2979 : * by calling the "finish_timestep_elem" methods of all passed element
2980 : * discretizations on one given subset.
2981 : *
2982 : * \param[in] vElemDisc element discretizations
2983 : * \param[in] dd DoF Distribution
2984 : * \param[in] si subset index
2985 : * \param[in] bNonRegularGrid flag to indicate if non regular grid is used
2986 : * \param[in] vSol current and previous solutions
2987 : */
2988 : template <typename TDomain, typename TAlgebra, typename TGlobAssembler>
2989 : template <typename TElem>
2990 0 : void DomainDiscretizationBase<TDomain, TAlgebra, TGlobAssembler>::
2991 : FinishTimestepElem(const std::vector<IElemDisc<domain_type>*>& vElemDisc,
2992 : ConstSmartPtr<DoFDistribution> dd,
2993 : int si, bool bNonRegularGrid,
2994 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol)
2995 : {
2996 : // check if only some elements are selected
2997 0 : if(m_spAssTuner->selected_elements_used())
2998 : {
2999 : std::vector<TElem*> vElem;
3000 0 : m_spAssTuner->collect_selected_elements(vElem, dd, si);
3001 :
3002 : // assembling is carried out only over those elements
3003 : // which are selected and in subset si
3004 : gass_type::template FinishTimestepElem<TElem>
3005 0 : (vElemDisc, m_spApproxSpace->domain(), dd, vElem.begin(), vElem.end(), si,
3006 : bNonRegularGrid, vSol, m_spAssTuner);
3007 0 : }
3008 : else
3009 : {
3010 : // general case: assembling over all elements in subset si
3011 : gass_type::template FinishTimestepElem<TElem>
3012 0 : (vElemDisc, m_spApproxSpace->domain(), dd,
3013 : dd->template begin<TElem>(si), dd->template end<TElem>(si), si,
3014 : bNonRegularGrid, vSol, m_spAssTuner);
3015 : }
3016 0 : }
3017 :
3018 : } // end namespace ug
3019 :
3020 : #endif /*__H__UG__LIB_DISC__SPATIAL_DISC__DOMAIN_DISC_IMPL__*/
|
