Multiscale Universal Interface  2.0
A Concurrent Framework for Coupling Heterogeneous Solvers
algo_aitken.h
Go to the documentation of this file.
1 /*****************************************************************************
2 * Multiscale Universal Interface Code Coupling Library *
3 * *
4 * Copyright (C) 2019 Y. H. Tang, S. Kudo, X. Bian, Z. Li, G. E. Karniadakis *
5 * W. Liu *
6 * *
7 * This software is jointly licensed under the Apache License, Version 2.0 *
8 * and the GNU General Public License version 3, you may use it according *
9 * to either. *
10 * *
11 * ** Apache License, version 2.0 ** *
12 * *
13 * Licensed under the Apache License, Version 2.0 (the "License"); *
14 * you may not use this file except in compliance with the License. *
15 * You may obtain a copy of the License at *
16 * *
17 * http://www.apache.org/licenses/LICENSE-2.0 *
18 * *
19 * Unless required by applicable law or agreed to in writing, software *
20 * distributed under the License is distributed on an "AS IS" BASIS, *
21 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. *
22 * See the License for the specific language governing permissions and *
23 * limitations under the License. *
24 * *
25 * ** GNU General Public License, version 3 ** *
26 * *
27 * This program is free software: you can redistribute it and/or modify *
28 * it under the terms of the GNU General Public License as published by *
29 * the Free Software Foundation, either version 3 of the License, or *
30 * (at your option) any later version. *
31 * *
32 * This program is distributed in the hope that it will be useful, *
33 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
34 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
35 * GNU General Public License for more details. *
36 * *
37 * You should have received a copy of the GNU General Public License *
38 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
39 *****************************************************************************/
40 
48 #ifndef MUI_COUPLING_ALGORITHM_AITKEN_H_
49 #define MUI_COUPLING_ALGORITHM_AITKEN_H_
50 
51 #include "../../general/util.h"
52 #include "../../config.h"
53 #include <mpi.h>
54 #include <iterator>
55 
56 namespace mui {
57 
58 template<typename CONFIG=default_config> class algo_aitken {
59 public:
60  using REAL = typename CONFIG::REAL;
61  using INT = typename CONFIG::INT;
62  using time_type = typename CONFIG::time_type;
63  using iterator_type = typename CONFIG::iterator_type;
64  using point_type = typename CONFIG::point_type;
65 
66  algo_aitken( REAL under_relaxation_factor = 1.0,
67  REAL under_relaxation_factor_max = 1.0,
68  MPI_Comm local_comm = MPI_COMM_NULL,
69  std::vector<std::pair<point_type, REAL>> pts_vlu_init =
70  std::vector<std::pair<point_type, REAL>>(),
71  REAL res_l2_norm_nm1 = 0.0):
72  init_under_relaxation_factor_(under_relaxation_factor),
73  under_relaxation_factor_max_(under_relaxation_factor_max),
74  local_mpi_comm_world_(local_comm) {
75  minimum_iterator_ = 1;
76  under_relaxation_factor_.insert(under_relaxation_factor_.begin(),
77  std::make_pair(std::make_pair(
78  std::numeric_limits<time_type>::lowest(),
79  (minimum_iterator_ -1)), init_under_relaxation_factor_
80  )
81  );
82 
83  if (!pts_vlu_init.empty()) {
84  pts_time_value_.insert(pts_time_value_.begin(),
85  std::make_pair(std::make_pair(
86  std::numeric_limits<time_type>::lowest(),
87  (minimum_iterator_ -1)),pts_vlu_init
88  )
89  );
90  }
91 
92  if(res_l2_norm_nm1 != 0.0){
93  residual_l2_norm_.insert(residual_l2_norm_.begin(),
94  std::make_pair(std::make_pair(
95  std::numeric_limits<time_type>::lowest(),
96  (minimum_iterator_ -1)),
97  (std::make_pair(static_cast<INT>(0), res_l2_norm_nm1)
98  )
99  )
100  );
101  }
102  }
103 
104  //- relaxation based on single time value
105  template<typename OTYPE>
106  OTYPE relaxation(std::pair<time_type,iterator_type> t, point_type focus, OTYPE filtered_value) const {
107 
108  OTYPE filtered_old_value = 0.0;
109 
110  if (pts_time_value_.empty()) {
111 
112  assert(pts_time_res_.empty());
113 
114  filtered_old_value = 0.0;
115 
116  std::vector<std::pair<point_type, REAL>> pts_value_temp{
117  std::make_pair(focus,calculate_relaxed_value(t,filtered_value,filtered_old_value))
118  };
119 
120  pts_time_value_.insert(pts_time_value_.begin(),std::make_pair(t,pts_value_temp));
121 
122  std::vector<std::pair<point_type, REAL>> pts_res_temp{
123  std::make_pair(focus,calculate_point_residual(t,filtered_value,filtered_old_value))
124  };
125 
126  pts_time_res_.insert(pts_time_res_.begin(),std::make_pair(t,pts_res_temp));
127 
128  return calculate_relaxed_value(t,filtered_value,filtered_old_value);
129 
130  } else {
131 
132  auto present_iter = std::find_if(pts_time_value_.begin(), pts_time_value_.end(),
133  [t](std::pair<std::pair<time_type,iterator_type>,std::vector<std::pair<point_type, REAL>>> b) {
134  return ((t.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
135  (t.second == b.first.second);});
136 
137  auto mi = minimum_iterator_;
138  auto previous_iter = std::find_if(pts_time_value_.begin(), pts_time_value_.end(),
139  [t, &mi](std::pair<std::pair<time_type,iterator_type>,std::vector<std::pair<point_type, REAL>>> b) {
140  return ((t.second == mi) ?
141  (b.first.first < t.first) ||
142  (((t.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
143  (b.first.second == (mi - 1))) :
144  ((t.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
145  (b.first.second < t.second));});
146 
147  auto present_res_iter = std::find_if(pts_time_res_.begin(), pts_time_res_.end(),
148  [t](std::pair<std::pair<time_type,iterator_type>,std::vector<std::pair<point_type, REAL>>> b) {
149  return ((t.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
150  (t.second == b.first.second);});
151 
152  mi = minimum_iterator_;
153  auto previous_res_iter = std::find_if(pts_time_res_.begin(), pts_time_res_.end(),
154  [t, &mi](std::pair<std::pair<time_type,iterator_type>,std::vector<std::pair<point_type, REAL>>> b) {
155  return ((t.second == mi) ?
156  (b.first.first < t.first) ||
157  (((t.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
158  (b.first.second == (mi - 1))) :
159  ((t.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
160  (b.first.second < t.second));});
161 
162  if ((present_iter == std::end(pts_time_value_)) &&
163  (previous_iter == std::end(pts_time_value_)) ) {
164 
165  assert((present_res_iter == std::end(pts_time_res_)) || pts_time_res_.empty());
166 
167  std::cerr << "Non-monotonic time marching does not (yet) supported for the Aitken coupling method! " << std::endl;
168 
169  } else if ((present_iter != std::end(pts_time_value_)) &&
170  (previous_iter == std::end(pts_time_value_)) ) {
171 
172  assert((((present_iter->first.first - present_res_iter->first.first) < std::numeric_limits<time_type>::epsilon()) &&
173  (present_iter->first.second == present_res_iter->first.second)) || pts_time_res_.empty());
174 
175  if (!pts_time_res_.empty()) {
176  assert(!residual_l2_norm_.empty());
177  }
178 
179  auto pts_relx_value_iter = std::find_if(present_iter->second.begin(),
180  present_iter->second.end(), [focus](std::pair<point_type, REAL> b) {
181  return normsq(focus - b.first) < std::numeric_limits<REAL>::epsilon();
182  });
183 
184  auto pts_relx_res_iter = std::find_if(present_res_iter->second.begin(),
185  present_res_iter->second.end(), [focus](std::pair<point_type, REAL> b) {
186  return normsq(focus - b.first) < std::numeric_limits<REAL>::epsilon();
187  });
188 
189  if ( pts_relx_value_iter == std::end(present_iter->second) ) {
190 
191  assert((pts_relx_res_iter == std::end(present_res_iter->second)) || pts_time_res_.empty());
192 
193  // Interpolate the relaxed value by N2_linear
194  REAL r2min_1st = std::numeric_limits<REAL>::max();
195  REAL r2min_2nd = std::numeric_limits<REAL>::max();
196  OTYPE value_1st = 0, value_2nd = 0;
197  for( size_t i = 0 ; i < present_iter->second.size() ; i++ ) {
198  REAL dr2 = normsq( focus - present_iter->second[i].first );
199  if ( dr2 < r2min_1st ) {
200  r2min_2nd = r2min_1st;
201  value_2nd = value_1st;
202  r2min_1st = dr2;
203  value_1st = present_iter->second[i].second ;
204  } else if ( dr2 < r2min_2nd ) {
205  r2min_2nd = dr2;
206  value_2nd = present_iter->second[i].second ;
207  }
208  }
209 
210  REAL r1 = std::sqrt( r2min_1st );
211  REAL r2 = std::sqrt( r2min_2nd );
212 
213  auto relaxed_value_temp = ( value_1st * r2 + value_2nd * r1 ) / ( r1 + r2 );
214 
215  present_iter->second.insert(present_iter->second.begin(),std::make_pair(focus, relaxed_value_temp));
216 
217  if (pts_time_res_.empty()) {
218 
219  std::vector<std::pair<point_type, REAL>> pts_res_temp{
220  std::make_pair(focus, (filtered_value-relaxed_value_temp))
221  };
222 
223  pts_time_res_.insert(pts_time_res_.begin(),std::make_pair(t,pts_res_temp));
224 
225  } else {
226  present_res_iter->second.insert(present_res_iter->second.begin(),std::make_pair(focus, (filtered_value-relaxed_value_temp)));
227  }
228 
229  return relaxed_value_temp;
230 
231  } else {
232 
233  assert((normsq(pts_relx_value_iter->first - pts_relx_res_iter->first) < std::numeric_limits<REAL>::epsilon()) || pts_time_res_.empty());
234 
235  return pts_relx_value_iter->second;
236 
237  }
238 
239  } else if ((present_iter == std::end(pts_time_value_)) &&
240  (previous_iter != std::end(pts_time_value_)) ) {
241 
242  assert((present_res_iter == std::end(pts_time_res_)) || pts_time_res_.empty());
243 
244  auto pts_relx_value_iter = std::find_if(previous_iter->second.begin(),
245  previous_iter->second.end(), [focus](std::pair<point_type, REAL> b) {
246  return normsq(focus - b.first) < std::numeric_limits<REAL>::epsilon();
247  });
248 
249  if ( pts_relx_value_iter == std::end(previous_iter->second) ) {
250 
251  // Interpolate the relaxed value by N2_linear
252  REAL r2min_1st = std::numeric_limits<REAL>::max();
253  REAL r2min_2nd = std::numeric_limits<REAL>::max();
254  OTYPE value_1st = 0, value_2nd = 0;
255  for( size_t i = 0 ; i < previous_iter->second.size() ; i++ ) {
256  REAL dr2 = normsq( focus - previous_iter->second[i].first );
257  if ( dr2 < r2min_1st ) {
258  r2min_2nd = r2min_1st;
259  value_2nd = value_1st;
260  r2min_1st = dr2;
261  value_1st = previous_iter->second[i].second ;
262  } else if ( dr2 < r2min_2nd ) {
263  r2min_2nd = dr2;
264  value_2nd = previous_iter->second[i].second ;
265  }
266  }
267 
268  REAL r1 = std::sqrt( r2min_1st );
269  REAL r2 = std::sqrt( r2min_2nd );
270 
271  filtered_old_value = ( value_1st * r2 + value_2nd * r1 ) / ( r1 + r2 );
272 
273  } else {
274 
275  filtered_old_value = pts_relx_value_iter->second;
276 
277  }
278 
279  std::vector<std::pair<point_type, REAL>> pts_value_temp{
280  std::make_pair(focus,calculate_relaxed_value(t,filtered_value,filtered_old_value))
281  };
282 
283  pts_time_value_.insert(pts_time_value_.begin(),std::make_pair(t, pts_value_temp));
284 
285  std::vector<std::pair<point_type, REAL>> pts_res_temp{
286  std::make_pair(focus,calculate_point_residual(t,filtered_value,filtered_old_value))
287  };
288 
289  pts_time_res_.insert(pts_time_res_.begin(),std::make_pair(t,pts_res_temp));
290 
291  present_iter = std::find_if(pts_time_value_.begin(), pts_time_value_.end(),
292  [t](std::pair<std::pair<time_type,iterator_type>,std::vector<std::pair<point_type, REAL>>> b) {
293  return ((t.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
294  (t.second == b.first.second);});
295 
296  auto mi = minimum_iterator_;
297  previous_iter = std::find_if(pts_time_value_.begin(), pts_time_value_.end(),
298  [t, &mi](std::pair<std::pair<time_type,iterator_type>,std::vector<std::pair<point_type, REAL>>> b) {
299  return ((t.second == mi) ?
300  (b.first.first < t.first) ||
301  (((t.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
302  (b.first.second == (mi - 1))) :
303  ((t.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
304  (b.first.second < t.second));});
305 
306  present_res_iter = std::find_if(pts_time_res_.begin(), pts_time_res_.end(),
307  [t](std::pair<std::pair<time_type,iterator_type>,std::vector<std::pair<point_type, REAL>>> b) {
308  return ((t.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
309  (t.second == b.first.second);});
310 
311  mi = minimum_iterator_;
312  previous_res_iter = std::find_if(pts_time_res_.begin(), pts_time_res_.end(),
313  [t, &mi](std::pair<std::pair<time_type,iterator_type>,std::vector<std::pair<point_type, REAL>>> b) {
314  return ((t.second == mi) ?
315  (b.first.first < t.first) ||
316  (((t.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
317  (b.first.second == (mi - 1))) :
318  ((t.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
319  (b.first.second < t.second));});
320 
321  auto pts_residual_l2_norm_iter = std::find_if(residual_l2_norm_.begin(),
322  residual_l2_norm_.end(), [previous_iter](std::pair<std::pair<time_type,iterator_type>,std::pair<INT, REAL>> b) {
323  return ((previous_iter->first.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
324  (previous_iter->first.second == b.first.second);});
325 
326  if(pts_residual_l2_norm_iter == std::end(residual_l2_norm_)) {
327 
328  if (previous_res_iter!=std::end(pts_time_res_)) {
329 
330  assert(((previous_iter->first.first - previous_res_iter->first.first) < std::numeric_limits<time_type>::epsilon()) &&
331  (previous_iter->first.second == previous_res_iter->first.second) );
332 
333  REAL local_residual_mag_sq_sum_temp = 0.0;
334  REAL residual_mag_sq_sum_temp = 0.0;
335 
336  for (auto & element_pair : previous_res_iter->second) {
337  local_residual_mag_sq_sum_temp += std::pow(element_pair.second, 2);
338  }
339 
340  if ( local_mpi_comm_world_ == MPI_COMM_NULL ) {
341  residual_mag_sq_sum_temp = local_residual_mag_sq_sum_temp;
342  } else {
343  MPI_Allreduce(&local_residual_mag_sq_sum_temp, &residual_mag_sq_sum_temp, 1, MPI_DOUBLE, MPI_SUM, local_mpi_comm_world_);
344  }
345 
346  if((residual_mag_sq_sum_temp != 0) || (!residual_l2_norm_.empty())){
347  residual_l2_norm_.insert(residual_l2_norm_.begin(),
348  std::make_pair(
349  previous_res_iter->first, (
350  std::make_pair(
351  static_cast<INT>(
352  previous_res_iter->second.size()
353  ), std::sqrt(residual_mag_sq_sum_temp)
354  )
355  )
356  )
357  );
358 
359  pts_residual_l2_norm_iter = std::find_if(residual_l2_norm_.begin(),
360  residual_l2_norm_.end(), [previous_iter](std::pair<std::pair<time_type,iterator_type>,std::pair<INT, REAL>> b) {
361  return ((previous_iter->first.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
362  (previous_iter->first.second == b.first.second);});
363  }
364  }
365 
366  } else {
367 
368  if(pts_residual_l2_norm_iter->second.first != 0) {
369 
370  auto pts_time_res_iter = std::find_if(pts_time_res_.begin(),
371  pts_time_res_.end(), [previous_res_iter](std::pair<std::pair<time_type,iterator_type>,std::vector<std::pair<point_type, REAL>>> b) {
372  return ((previous_res_iter->first.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
373  (previous_res_iter->first.second == b.first.second);});
374 
375  assert(pts_time_res_iter != std::end(pts_time_res_) );
376 
377  if(pts_time_res_iter->second.size() != static_cast<size_t>(pts_residual_l2_norm_iter->second.first)) {
378 
379  REAL local_residual_mag_sq_sum_temp = 0.0;
380  REAL residual_mag_sq_sum_temp = 0.0;
381 
382  for (auto & element_pair : pts_time_res_iter->second) {
383  local_residual_mag_sq_sum_temp += std::pow(element_pair.second, 2);
384  }
385 
386  if ( local_mpi_comm_world_ == MPI_COMM_NULL ) {
387  residual_mag_sq_sum_temp = local_residual_mag_sq_sum_temp;
388  } else {
389  MPI_Allreduce(&local_residual_mag_sq_sum_temp, &residual_mag_sq_sum_temp, 1, MPI_DOUBLE, MPI_SUM, local_mpi_comm_world_);
390  }
391 
392  pts_residual_l2_norm_iter->second.second = std::sqrt(residual_mag_sq_sum_temp);
393  }
394  }
395  }
396 
397  return calculate_relaxed_value(t,filtered_value,filtered_old_value);
398 
399  } else {
400 
401  assert((((present_iter->first.first - present_res_iter->first.first) < std::numeric_limits<time_type>::epsilon()) &&
402  (present_iter->first.second == present_res_iter->first.second)) || pts_time_res_.empty());
403 
404  auto pts_relx_value_iter = std::find_if(previous_iter->second.begin(),
405  previous_iter->second.end(), [focus](std::pair<point_type, REAL> b) {
406  return normsq(focus - b.first) < std::numeric_limits<REAL>::epsilon();
407  });
408 
409  if ( pts_relx_value_iter == std::end(previous_iter->second) ) {
410 
411  // Interpolate the relaxed value by N2_linear
412  REAL r2min_1st = std::numeric_limits<REAL>::max();
413  REAL r2min_2nd = std::numeric_limits<REAL>::max();
414  OTYPE value_1st = 0, value_2nd = 0;
415  for( size_t i = 0 ; i < previous_iter->second.size() ; i++ ) {
416  REAL dr2 = normsq( focus - previous_iter->second[i].first );
417  if ( dr2 < r2min_1st ) {
418  r2min_2nd = r2min_1st;
419  value_2nd = value_1st;
420  r2min_1st = dr2;
421  value_1st = previous_iter->second[i].second ;
422  } else if ( dr2 < r2min_2nd ) {
423  r2min_2nd = dr2;
424  value_2nd = previous_iter->second[i].second ;
425  }
426  }
427 
428  REAL r1 = std::sqrt( r2min_1st );
429  REAL r2 = std::sqrt( r2min_2nd );
430 
431  filtered_old_value = ( value_1st * r2 + value_2nd * r1 ) / ( r1 + r2 );
432 
433  } else {
434 
435  filtered_old_value = pts_relx_value_iter->second;
436 
437  }
438 
439  auto pts_present_relx_value_iter = std::find_if(present_iter->second.begin(),
440  present_iter->second.end(), [focus](std::pair<point_type, REAL> b) {
441  return normsq(focus - b.first) < std::numeric_limits<REAL>::epsilon();
442  });
443 
444  auto pts_relx_res_iter = std::find_if(present_res_iter->second.begin(),
445  present_res_iter->second.end(), [focus](std::pair<point_type, REAL> b) {
446  return normsq(focus - b.first) < std::numeric_limits<REAL>::epsilon();
447  });
448 
449  if ( pts_present_relx_value_iter == std::end(present_iter->second) ) {
450 
451  assert((pts_relx_res_iter == std::end(present_res_iter->second)) || pts_time_res_.empty());
452 
453  present_iter->second.insert(present_iter->second.begin(),
454  std::make_pair(focus,calculate_relaxed_value(
455  t,filtered_value,filtered_old_value)
456  )
457  );
458 
459  if (pts_time_res_.empty()) {
460 
461  std::vector<std::pair<point_type, REAL>> pts_res_temp{
462  std::make_pair(focus,calculate_point_residual(
463  t,filtered_value,filtered_old_value))
464  };
465 
466  pts_time_res_.insert(pts_time_res_.begin(),std::make_pair(t,pts_res_temp));
467 
468  } else {
469 
470  present_res_iter->second.insert(present_res_iter->second.begin(),
471  std::make_pair(focus,calculate_point_residual(
472  t,filtered_value,filtered_old_value)
473  )
474  );
475  }
476  } else {
477 
478  assert((normsq(pts_present_relx_value_iter->first - pts_relx_res_iter->first) < std::numeric_limits<REAL>::epsilon()) || pts_time_res_.empty());
479 
480  pts_present_relx_value_iter->second = calculate_relaxed_value(t,filtered_value,filtered_old_value);
481 
482  if (pts_time_res_.empty()) {
483 
484  std::vector<std::pair<point_type, REAL>> pts_res_temp{
485  std::make_pair(focus,calculate_point_residual(
486  t,filtered_value,filtered_old_value))
487  };
488 
489  pts_time_res_.insert(pts_time_res_.begin(),std::make_pair(t,pts_res_temp));
490 
491  } else {
492  pts_relx_res_iter->second = calculate_point_residual(t,filtered_value,filtered_old_value);
493  }
494  }
495 
496  return calculate_relaxed_value(t,filtered_value,filtered_old_value);
497  }
498  return calculate_relaxed_value(t,filtered_value,filtered_old_value);
499  }
500  }
501 
502  REAL get_under_relaxation_factor (time_type t_single) {
503 
504  std::pair<time_type,iterator_type> t = std::make_pair(std::numeric_limits<time_type>::lowest(),t_single);
505 
506  update_under_relaxation_factor(t);
507 
508  auto under_relaxation_present_iter = std::find_if(under_relaxation_factor_.begin(),
509  under_relaxation_factor_.end(), [t](std::pair<std::pair<time_type,iterator_type>, REAL> b) {
510  return ((t.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
511  (t.second == b.first.second);});
512 
513  assert(under_relaxation_present_iter != std::end(under_relaxation_factor_) );
514 
515  return under_relaxation_present_iter->second;
516 
517  }
518 
519  REAL get_under_relaxation_factor (time_type t, iterator_type it) {
520 
521  std::pair<time_type,iterator_type> time = std::make_pair(t,it);
522 
523  update_under_relaxation_factor(time);
524 
525  auto under_relaxation_present_iter = std::find_if(under_relaxation_factor_.begin(),
526  under_relaxation_factor_.end(), [time](std::pair<std::pair<time_type,iterator_type>, REAL> b) {
527  return ((time.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
528  (time.second == b.first.second);});
529 
530  assert(under_relaxation_present_iter != std::end(under_relaxation_factor_) );
531 
532  return under_relaxation_present_iter->second;
533 
534  }
535 
536  REAL get_residual_L2_Norm (time_type t_single) {
537 
538  std::pair<time_type,iterator_type> t = std::make_pair(std::numeric_limits<time_type>::lowest(),t_single);
539 
540  update_under_relaxation_factor(t);
541 
542  auto mi = minimum_iterator_;
543  auto res_l2_norm_nm1_iter = std::find_if(residual_l2_norm_.begin(),
544  residual_l2_norm_.end(), [t, &mi](std::pair<std::pair<time_type,iterator_type>,std::pair<INT, REAL>> b) {
545  return ((t.second == mi) ?
546  (b.first.first < t.first) ||
547  (((t.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
548  (b.first.second == (mi - 1))) :
549  ((t.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
550  (b.first.second < t.second));});
551 
552  if(res_l2_norm_nm1_iter != std::end(residual_l2_norm_) ) {
553  return res_l2_norm_nm1_iter->second.second;
554  } else {
555  return 0.0;
556  }
557  }
558 
559  REAL get_residual_L2_Norm (time_type t, iterator_type it) {
560 
561  std::pair<time_type,iterator_type> time = std::make_pair(t,it);
562 
563  update_under_relaxation_factor(time);
564 
565  auto mi = minimum_iterator_;
566  auto res_l2_norm_nm1_iter = std::find_if(residual_l2_norm_.begin(),
567  residual_l2_norm_.end(), [time, &mi](std::pair<std::pair<time_type,iterator_type>,std::pair<INT, REAL>> b) {
568  return ((time.second == mi) ?
569  (b.first.first < time.first) ||
570  (((time.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
571  (b.first.second == (mi - 1))) :
572  ((time.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
573  (b.first.second < time.second));});
574 
575  if(res_l2_norm_nm1_iter != std::end(residual_l2_norm_) ) {
576  return res_l2_norm_nm1_iter->second.second;
577  } else {
578  return 0.0;
579  }
580  }
581 
582 private:
583  template<typename OTYPE>
584  OTYPE calculate_relaxed_value(std::pair<time_type,iterator_type> t, OTYPE filtered_value, OTYPE filtered_old_value) const {
585 
586  update_under_relaxation_factor(t);
587 
588  auto under_relaxation_present_iter = std::find_if(under_relaxation_factor_.begin(),
589  under_relaxation_factor_.end(), [t](std::pair<std::pair<time_type,iterator_type>, REAL> b) {
590  return ((t.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
591  (t.second == b.first.second);});
592 
593  assert(under_relaxation_present_iter != std::end(under_relaxation_factor_) );
594 
595  return (under_relaxation_present_iter->second * filtered_value) +
596  ((1 - under_relaxation_present_iter->second) * filtered_old_value);
597 
598  }
599 
600  void update_under_relaxation_factor(std::pair<time_type,iterator_type> t) const {
601 
602  auto under_relaxation_present_iter = std::find_if(under_relaxation_factor_.begin(),
603  under_relaxation_factor_.end(), [t](std::pair<std::pair<time_type,iterator_type>, REAL> b) {
604  return ((t.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
605  (t.second == b.first.second);});
606 
607  auto mi = minimum_iterator_;
608  auto under_relaxation_prev_iter = std::find_if(under_relaxation_factor_.begin(),
609  under_relaxation_factor_.end(), [t, &mi](std::pair<std::pair<time_type,iterator_type>, REAL> b) {
610  return ((t.second == mi) ?
611  (b.first.first < t.first) ||
612  (((t.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
613  (b.first.second == (mi - 1))) :
614  ((t.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
615  (b.first.second < t.second));});
616 
617  mi = minimum_iterator_;
618  auto res_l2_norm_nm1_iter = std::find_if(residual_l2_norm_.begin(),
619  residual_l2_norm_.end(), [t, &mi](std::pair<std::pair<time_type,iterator_type>,std::pair<INT, REAL>> b) {
620  return ((t.second == mi) ?
621  (b.first.first < t.first) ||
622  (((t.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
623  (b.first.second == (mi - 1))) :
624  ((t.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
625  (b.first.second < t.second));});
626 
627  mi = minimum_iterator_;
628  auto res_l2_norm_nm2_iter = std::find_if(residual_l2_norm_.begin(),
629  residual_l2_norm_.end(), [res_l2_norm_nm1_iter, &mi](std::pair<std::pair<time_type,iterator_type>,std::pair<INT, REAL>> b) {
630  return ((res_l2_norm_nm1_iter->first.second == mi) ?
631  (b.first.first < res_l2_norm_nm1_iter->first.first) ||
632  (((res_l2_norm_nm1_iter->first.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
633  (b.first.second == (mi - 1))) :
634  ((res_l2_norm_nm1_iter->first.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
635  (b.first.second < res_l2_norm_nm1_iter->first.second));});
636 
637 
638  if (t.second <= minimum_iterator_) {
639 
640  under_relaxation_factor_.insert(under_relaxation_factor_.begin(),std::make_pair(t, init_under_relaxation_factor_));
641 
642  under_relaxation_present_iter = std::find_if(under_relaxation_factor_.begin(),
643  under_relaxation_factor_.end(), [t](std::pair<std::pair<time_type,iterator_type>, REAL> b) {
644  return ((t.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
645  (t.second == b.first.second);});
646 
647  auto mi = minimum_iterator_;
648  under_relaxation_prev_iter = std::find_if(under_relaxation_factor_.begin(),
649  under_relaxation_factor_.end(), [t, &mi](std::pair<std::pair<time_type,iterator_type>, REAL> b) {
650  return ((t.second == mi) ?
651  (b.first.first < t.first) ||
652  (((t.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
653  (b.first.second == (mi - 1))) :
654  ((t.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
655  (b.first.second < t.second));});
656 
657  } else {
658 
659  if (under_relaxation_present_iter==std::end(under_relaxation_factor_)) {
660 
661  if(res_l2_norm_nm2_iter == std::end(residual_l2_norm_) ) {
662 
663  under_relaxation_factor_.insert(under_relaxation_factor_.begin(),std::make_pair(t, init_under_relaxation_factor_));
664 
665  under_relaxation_present_iter = std::find_if(under_relaxation_factor_.begin(),
666  under_relaxation_factor_.end(), [t](std::pair<std::pair<time_type,iterator_type>, REAL> b) {
667  return ((t.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
668  (t.second == b.first.second);});
669 
670  auto mi = minimum_iterator_;
671  under_relaxation_prev_iter = std::find_if(under_relaxation_factor_.begin(),
672  under_relaxation_factor_.end(), [t, &mi](std::pair<std::pair<time_type,iterator_type>, REAL> b) {
673  return ((t.second == mi) ?
674  (b.first.first < t.first) ||
675  (((t.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
676  (b.first.second == (mi - 1))) :
677  ((t.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
678  (b.first.second < t.second));});
679 
680 
681  } else {
682 
683  if(res_l2_norm_nm2_iter->second.first != 0 ) {
684  auto pts_time_res_nm2_iter = std::find_if(pts_time_res_.begin(),
685  pts_time_res_.end(), [res_l2_norm_nm2_iter](std::pair<std::pair<time_type,iterator_type>,std::vector<std::pair<point_type, REAL>>> b) {
686  return ((res_l2_norm_nm2_iter->first.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
687  (res_l2_norm_nm2_iter->first.second == b.first.second);});
688 
689  assert(pts_time_res_nm2_iter != std::end(pts_time_res_) );
690 
691  if(pts_time_res_nm2_iter->second.size() != static_cast<size_t>(res_l2_norm_nm2_iter->second.first)) {
692 
693  REAL local_residual_mag_sq_sum_temp = 0.0;
694  REAL residual_mag_sq_sum_temp = 0.0;
695 
696  for (auto & element_pair : pts_time_res_nm2_iter->second) {
697  local_residual_mag_sq_sum_temp += std::pow(element_pair.second, 2);
698  }
699 
700  if ( local_mpi_comm_world_ == MPI_COMM_NULL ) {
701  residual_mag_sq_sum_temp = local_residual_mag_sq_sum_temp;
702  } else {
703  MPI_Allreduce(&local_residual_mag_sq_sum_temp, &residual_mag_sq_sum_temp, 1, MPI_DOUBLE, MPI_SUM, local_mpi_comm_world_);
704  }
705 
706  res_l2_norm_nm2_iter->second.second = std::sqrt(residual_mag_sq_sum_temp);
707 
708  update_under_relaxation_factor(res_l2_norm_nm2_iter->first);
709  }
710 
711  auto pts_time_res_nm1_iter = std::find_if(pts_time_res_.begin(),
712  pts_time_res_.end(), [res_l2_norm_nm1_iter](std::pair<std::pair<time_type,iterator_type>,std::vector<std::pair<point_type, REAL>>> b) {
713  return ((res_l2_norm_nm1_iter->first.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
714  (res_l2_norm_nm1_iter->first.second == b.first.second);});
715 
716  if(pts_time_res_nm1_iter->second.size() != static_cast<size_t>(res_l2_norm_nm1_iter->second.first)) {
717 
718  REAL local_residual_mag_sq_sum_temp = 0.0;
719  REAL residual_mag_sq_sum_temp = 0.0;
720 
721  for (auto & element_pair : pts_time_res_nm1_iter->second) {
722  local_residual_mag_sq_sum_temp += std::pow(element_pair.second, 2);
723  }
724 
725  if ( local_mpi_comm_world_ == MPI_COMM_NULL ) {
726  residual_mag_sq_sum_temp = local_residual_mag_sq_sum_temp;
727  } else {
728  MPI_Allreduce(&local_residual_mag_sq_sum_temp, &residual_mag_sq_sum_temp, 1, MPI_DOUBLE, MPI_SUM, local_mpi_comm_world_);
729  }
730 
731  res_l2_norm_nm1_iter->second.second = std::sqrt(residual_mag_sq_sum_temp);
732 
733  update_under_relaxation_factor(res_l2_norm_nm1_iter->first);
734  }
735  }
736 
737  double nominator = res_l2_norm_nm2_iter->second.second;
738  double denominator = res_l2_norm_nm1_iter->second.second - res_l2_norm_nm2_iter->second.second;
739 
740  if (denominator != 0.0 ) {
741 
742  if (under_relaxation_prev_iter==std::end(under_relaxation_factor_)) {
743  under_relaxation_factor_.insert(under_relaxation_factor_.begin(),
744  std::make_pair(
745  t, calculate_aitken_constraint_pnz_control(init_under_relaxation_factor_)
746  )
747  );
748 
749  under_relaxation_present_iter = std::find_if(under_relaxation_factor_.begin(),
750  under_relaxation_factor_.end(), [t](std::pair<std::pair<time_type,iterator_type>, REAL> b) {
751  return ((t.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
752  (t.second == b.first.second);
753  });
754 
755  auto mi = minimum_iterator_;
756  under_relaxation_prev_iter = std::find_if(under_relaxation_factor_.begin(),
757  under_relaxation_factor_.end(), [t, &mi](std::pair<std::pair<time_type,iterator_type>, REAL> b) {
758  return ((t.second == mi) ?
759  (b.first.first < t.first) ||
760  (((t.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
761  (b.first.second == (mi - 1))) :
762  ((t.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
763  (b.first.second < t.second));});
764  } else {
765  under_relaxation_factor_.insert(under_relaxation_factor_.begin(),
766  std::make_pair(
767  t, calculate_aitken_constraint_pnz_control(
768  -under_relaxation_prev_iter->second * (nominator/denominator)
769  )
770  )
771  );
772 
773  under_relaxation_present_iter = std::find_if(under_relaxation_factor_.begin(),
774  under_relaxation_factor_.end(), [t](std::pair<std::pair<time_type,iterator_type>, REAL> b) {
775  return ((t.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
776  (t.second == b.first.second);});
777 
778  auto mi = minimum_iterator_;
779  under_relaxation_prev_iter = std::find_if(under_relaxation_factor_.begin(),
780  under_relaxation_factor_.end(), [t, &mi](std::pair<std::pair<time_type,iterator_type>, REAL> b) {
781  return ((t.second == mi) ?
782  (b.first.first < t.first) ||
783  (((t.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
784  (b.first.second == (mi - 1))) :
785  ((t.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
786  (b.first.second < t.second));});
787  }
788 
789  } else {
790  under_relaxation_factor_.insert(under_relaxation_factor_.begin(),std::make_pair(t, init_under_relaxation_factor_));
791 
792  under_relaxation_present_iter = std::find_if(under_relaxation_factor_.begin(),
793  under_relaxation_factor_.end(), [t](std::pair<std::pair<time_type,iterator_type>, REAL> b) {
794  return ((t.first - b.first.first) < std::numeric_limits<REAL>::epsilon()) &&
795  (t.second == b.first.second);});
796 
797  auto mi = minimum_iterator_;
798  under_relaxation_prev_iter = std::find_if(under_relaxation_factor_.begin(),
799  under_relaxation_factor_.end(), [t, &mi](std::pair<std::pair<time_type,iterator_type>, REAL> b) {
800  return ((t.second == mi) ?
801  (b.first.first < t.first) ||
802  (((t.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
803  (b.first.second == (mi - 1))) :
804  ((t.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
805  (b.first.second < t.second));});
806  }
807  }
808  } else {
809 
810  assert(under_relaxation_present_iter!=std::end(under_relaxation_factor_));
811 
812  if(res_l2_norm_nm2_iter == std::end(residual_l2_norm_) ) {
813  if(under_relaxation_present_iter->second != init_under_relaxation_factor_) {
814  std::cout << "change under Relax Factor to its initial value." << std::endl;
815  under_relaxation_present_iter->second = init_under_relaxation_factor_;
816  }
817  } else {
818 
819  if(res_l2_norm_nm2_iter->second.first != 0 ) {
820  auto pts_time_res_nm2_iter = std::find_if(pts_time_res_.begin(),
821  pts_time_res_.end(), [res_l2_norm_nm2_iter](std::pair<std::pair<time_type,iterator_type>,std::vector<std::pair<point_type, REAL>>> b) {
822  return ((res_l2_norm_nm2_iter->first.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
823  (res_l2_norm_nm2_iter->first.second == b.first.second);});
824 
825  assert(pts_time_res_nm2_iter != std::end(pts_time_res_) );
826 
827  if(pts_time_res_nm2_iter->second.size() != static_cast<size_t>(res_l2_norm_nm2_iter->second.first)) {
828 
829  REAL local_residual_mag_sq_sum_temp = 0.0;
830  REAL residual_mag_sq_sum_temp = 0.0;
831 
832  for (auto & element_pair : pts_time_res_nm2_iter->second) {
833  local_residual_mag_sq_sum_temp += std::pow(element_pair.second, 2);
834  }
835 
836  if ( local_mpi_comm_world_ == MPI_COMM_NULL ) {
837  residual_mag_sq_sum_temp = local_residual_mag_sq_sum_temp;
838  } else {
839  MPI_Allreduce(&local_residual_mag_sq_sum_temp, &residual_mag_sq_sum_temp, 1, MPI_DOUBLE, MPI_SUM, local_mpi_comm_world_);
840  }
841 
842  res_l2_norm_nm2_iter->second.second = std::sqrt(residual_mag_sq_sum_temp);
843 
844  }
845 
846  auto pts_time_res_nm1_iter = std::find_if(pts_time_res_.begin(),
847  pts_time_res_.end(), [res_l2_norm_nm1_iter](std::pair<std::pair<time_type,iterator_type>,std::vector<std::pair<point_type, REAL>>> b) {
848  return (((res_l2_norm_nm1_iter->first.first - b.first.first) < std::numeric_limits<time_type>::epsilon()) &&
849  (res_l2_norm_nm1_iter->first.second == b.first.second));});
850 
851  if(pts_time_res_nm1_iter->second.size() != static_cast<size_t>(res_l2_norm_nm1_iter->second.first)) {
852 
853  REAL local_residual_mag_sq_sum_temp = 0.0;
854  REAL residual_mag_sq_sum_temp = 0.0;
855 
856  for (auto & element_pair : pts_time_res_nm1_iter->second) {
857  local_residual_mag_sq_sum_temp += std::pow(element_pair.second, 2);
858  }
859 
860  if ( local_mpi_comm_world_ == MPI_COMM_NULL ) {
861  residual_mag_sq_sum_temp = local_residual_mag_sq_sum_temp;
862  } else {
863  MPI_Allreduce(&local_residual_mag_sq_sum_temp, &residual_mag_sq_sum_temp, 1, MPI_DOUBLE, MPI_SUM, local_mpi_comm_world_);
864  }
865 
866  res_l2_norm_nm1_iter->second.second = std::sqrt(residual_mag_sq_sum_temp);
867 
868  update_under_relaxation_factor(res_l2_norm_nm1_iter->first);
869  }
870  }
871 
872  double nominator = res_l2_norm_nm2_iter->second.second;
873  double denominator = res_l2_norm_nm1_iter->second.second - res_l2_norm_nm2_iter->second.second;
874 
875  if (denominator != 0.0 ) {
876 
877  if (under_relaxation_prev_iter==std::end(under_relaxation_factor_)) {
878  if(under_relaxation_present_iter->second != calculate_aitken_constraint_pnz_control(init_under_relaxation_factor_)
879  ) {
880  std::cout << "Update under Relx Factor. Position 01." << std::endl;
881  under_relaxation_present_iter->second = calculate_aitken_constraint_pnz_control(init_under_relaxation_factor_);
882  }
883  } else {
884  if(under_relaxation_present_iter->second != calculate_aitken_constraint_pnz_control(
885  -under_relaxation_prev_iter->second * (nominator/denominator))
886  ) {
887  std::cout << "Update under Relx Factor. Position 02." << std::endl;
888  under_relaxation_present_iter->second = calculate_aitken_constraint_pnz_control(
889  -under_relaxation_prev_iter->second * (nominator/denominator)
890  );
891  }
892  }
893  } else {
894  if(under_relaxation_present_iter->second != init_under_relaxation_factor_) {
895  std::cout << "Update under Relx Factor. Position 03." << std::endl;
896  under_relaxation_present_iter->second = init_under_relaxation_factor_;
897  }
898  }
899  }
900  }
901  }
902  }
903 
904  REAL calculate_aitken_constraint(REAL under_relaxation_factor) {
905 
906  return (sign(under_relaxation_factor) * std::min(std::abs(under_relaxation_factor),under_relaxation_factor_max_));
907 
908  }
909 
910  REAL calculate_aitken_constraint_pn_control(REAL under_relaxation_factor) {
911 
912  return (std::min(std::abs(under_relaxation_factor),under_relaxation_factor_max_));
913  }
914 
915  REAL calculate_aitken_constraint_pnz_control(REAL under_relaxation_factor) const {
916 
917  return ((std::min(std::abs(under_relaxation_factor),under_relaxation_factor_max_)) < init_under_relaxation_factor_) ? init_under_relaxation_factor_ : (std::min(std::abs(under_relaxation_factor),under_relaxation_factor_max_));
918  }
919 
920  template<typename OTYPE>
921  OTYPE calculate_point_residual(std::pair<time_type,iterator_type> t, OTYPE filtered_value, OTYPE filtered_old_value) const {
922 
923  return (filtered_value - calculate_relaxed_value(t, filtered_value, filtered_old_value));
924 
925  }
926 
927  INT sign(REAL value)
928  {
929  return (value < 0) ? -1 : ((value > 0) ? 1 : 0);
930  }
931 
932 protected:
933  const REAL init_under_relaxation_factor_;
934 
935  const REAL under_relaxation_factor_max_;
936 
937  MPI_Comm local_mpi_comm_world_;
938 
939  iterator_type minimum_iterator_;
940 
941  mutable std::vector<std::pair<std::pair<time_type,iterator_type>, REAL>> under_relaxation_factor_;
942 
943  mutable std::vector<std::pair<std::pair<time_type,iterator_type>,std::pair<INT, REAL>>> residual_l2_norm_;
944 
945  mutable std::vector<std::pair<std::pair<time_type,iterator_type>,std::vector<std::pair<point_type, REAL>>>> pts_time_value_;
946 
947  mutable std::vector<std::pair<std::pair<time_type,iterator_type>,std::vector<std::pair<point_type, REAL>>>> pts_time_res_;
948 
949 };
950 
951 }
952 
953 #endif /* MUI_COUPLING_ALGORITHM_AITKEN_H_ */
mui::algo_aitken
Definition: algo_aitken.h:58
mui::algo_aitken::get_residual_L2_Norm
REAL get_residual_L2_Norm(time_type t_single)
Definition: algo_aitken.h:536
mui::algo_aitken::minimum_iterator_
iterator_type minimum_iterator_
Definition: algo_aitken.h:939
mui::algo_aitken::pts_time_value_
std::vector< std::pair< std::pair< time_type, iterator_type >, std::vector< std::pair< point_type, REAL > > > > pts_time_value_
Definition: algo_aitken.h:945
mui::algo_aitken::get_residual_L2_Norm
REAL get_residual_L2_Norm(time_type t, iterator_type it)
Definition: algo_aitken.h:559
mui::algo_aitken::init_under_relaxation_factor_
const REAL init_under_relaxation_factor_
Definition: algo_aitken.h:933
mui::algo_aitken::under_relaxation_factor_max_
const REAL under_relaxation_factor_max_
Definition: algo_aitken.h:935
mui::algo_aitken::REAL
typename CONFIG::REAL REAL
Definition: algo_aitken.h:60
mui::algo_aitken::get_under_relaxation_factor
REAL get_under_relaxation_factor(time_type t, iterator_type it)
Definition: algo_aitken.h:519
mui::algo_aitken::under_relaxation_factor_
std::vector< std::pair< std::pair< time_type, iterator_type >, REAL > > under_relaxation_factor_
Definition: algo_aitken.h:941
mui::algo_aitken::relaxation
OTYPE relaxation(std::pair< time_type, iterator_type > t, point_type focus, OTYPE filtered_value) const
Definition: algo_aitken.h:106
mui::algo_aitken::INT
typename CONFIG::INT INT
Definition: algo_aitken.h:61
mui::algo_aitken::point_type
typename CONFIG::point_type point_type
Definition: algo_aitken.h:64
mui::algo_aitken::get_under_relaxation_factor
REAL get_under_relaxation_factor(time_type t_single)
Definition: algo_aitken.h:502
mui::algo_aitken::algo_aitken
algo_aitken(REAL under_relaxation_factor=1.0, REAL under_relaxation_factor_max=1.0, MPI_Comm local_comm=MPI_COMM_NULL, std::vector< std::pair< point_type, REAL >> pts_vlu_init=std::vector< std::pair< point_type, REAL >>(), REAL res_l2_norm_nm1=0.0)
Definition: algo_aitken.h:66
mui::algo_aitken::local_mpi_comm_world_
MPI_Comm local_mpi_comm_world_
Definition: algo_aitken.h:937
mui::algo_aitken::time_type
typename CONFIG::time_type time_type
Definition: algo_aitken.h:62
mui::algo_aitken::residual_l2_norm_
std::vector< std::pair< std::pair< time_type, iterator_type >, std::pair< INT, REAL > > > residual_l2_norm_
Definition: algo_aitken.h:943
mui::algo_aitken::pts_time_res_
std::vector< std::pair< std::pair< time_type, iterator_type >, std::vector< std::pair< point_type, REAL > > > > pts_time_res_
Definition: algo_aitken.h:947
mui::algo_aitken::iterator_type
typename CONFIG::iterator_type iterator_type
Definition: algo_aitken.h:63
mui::dim::time
dim< 0, 0, 1, 0, 0, 0, 0, 0 > time
Definition: dim.h:207
mui::dim::min
u u u u u u min
Definition: dim.h:289
mui
Definition: comm.h:54
mui::max
SCALAR max(vexpr< E, SCALAR, D > const &u)
Definition: point.h:350
mui::normsq
SCALAR normsq(vexpr< E, SCALAR, D > const &u)
Definition: point.h:380