Hermes2D  3.0
filter.cpp
1 // This file is part of Hermes2D.
2 //
3 // Hermes2D is free software: you can redistribute it and/or modify
4 // it under the terms of the GNU General Public License as published by
5 // the Free Software Foundation, either version 2 of the License, or
6 // (at your option) any later version.
7 //
8 // Hermes2D is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 // GNU General Public License for more details.
12 //
13 // You should have received a copy of the GNU General Public License
14 // along with Hermes2D. If not, see <http://www.gnu.org/licenses/>.
15 
16 #include "filter.h"
17 #include "forms.h"
18 #include "quad.h"
19 #include "refmap.h"
20 #include "traverse.h"
21 
22 namespace Hermes
23 {
24  namespace Hermes2D
25  {
26  template<typename Scalar>
27  Filter<Scalar>::Filter()
28  {
29  }
30 
31  template<typename Scalar>
32  Filter<Scalar>::Filter(std::vector<MeshFunctionSharedPtr<Scalar> > solutions) : MeshFunction<Scalar>(), solutions(solutions)
33  {
34  this->init();
35  }
36 
37  template<typename Scalar>
39  {
40  // construct the union mesh, if necessary
41  std::vector<MeshSharedPtr> meshes;
42  for (int i = 0; i < this->solutions.size(); i++)
43  meshes.push_back(this->solutions[i]->get_mesh());
44  this->mesh = meshes[0];
45 
46  for (int i = 0; i < this->solutions.size(); i++)
47  this->solutions_sub_idx.push_back(0);
48 
49  Solution<Scalar>* sln = dynamic_cast<Solution<Scalar>*>(this->solutions[0].get());
50  if (sln == nullptr)
51  this->space_type = HERMES_INVALID_SPACE;
52  else
53  this->space_type = sln->get_space_type();
54 
55  unimesh = false;
56 
57  for (int i = 1; i < this->solutions.size(); i++)
58  {
59  if (meshes[i] == nullptr)
60  {
61  this->warn("You may be initializing a Filter with Solution that is missing a Mesh.");
62  throw Hermes::Exceptions::Exception("this->meshes[%d] is nullptr in Filter<Scalar>::init().", i);
63  }
64  if (meshes[i]->get_seq() != this->mesh->get_seq())
65  {
66  unimesh = true;
67  break;
68  }
69 
70  sln = dynamic_cast<Solution<Scalar>*>(this->solutions[i].get());
71  if (sln == nullptr || sln->get_space_type() != this->space_type)
72  this->space_type = HERMES_INVALID_SPACE;
73  }
74 
75  if (unimesh)
76  {
77  this->mesh = MeshSharedPtr(new Mesh);
78  this->unidata = Traverse::construct_union_mesh(this->solutions.size(), &meshes[0], this->mesh);
79  }
80 
81  // misc init
82  this->num_components = 1;
83  this->order = 0;
84 
85  set_quad_2d(&g_quad_2d_std);
86  }
87 
88  template<typename Scalar>
90  {
91  this->free();
92  }
93 
94  template<typename Scalar>
96  {
98  for (int i = 0; i < this->solutions.size(); i++)
99  // nodup
100  this->solutions[i]->set_quad_2d(quad_2d);
101  }
102 
103  template<typename Scalar>
105  {
107  if (!unimesh)
108  {
109  for (int i = 0; i < this->solutions.size(); i++)
110  // nodup
111  this->solutions[i]->set_active_element(e);
112 
113  for (int i = 0; i < this->solutions.size(); i++)
114  solutions_sub_idx[i] = 0;
115  }
116  else
117  {
118  for (int i = 0; i < this->solutions.size(); i++)
119  {
120  this->solutions[i]->set_active_element(unidata[i][e->id].e);
121  this->solutions[i]->set_transform(unidata[i][e->id].idx);
122  solutions_sub_idx[i] = this->solutions[i]->get_transform();
123  }
124  }
125 
126  // fixme
127  this->order = 20;
128  }
129 
130  template<typename Scalar>
132  {
133  if (unimesh)
134  {
135  for (int i = 0; i < this->solutions.size(); i++)
136  free_with_check(unidata[i]);
137  free_with_check(unidata);
138  }
139  }
140 
141  template<typename Scalar>
143  {
144  free();
145  init();
146  }
147 
148  template<typename Scalar>
150  {
152  for (int i = 0; i < this->solutions.size(); i++)
153  {
154  // solutions_sub_idx[i] contains the value sln[i]->sub_idx, which the Filter thinks
155  // the solution has, or at least had last time. If the filter graph is
156  // cyclic, it could happen that some solutions would get pushed the transform
157  // more than once. This mechanism prevents it. If the filter sees that the
158  // solution already has a different sub_idx than it thinks it should have,
159  // it assumes someone else has already pushed the correct transform. This
160  // is actually the case for cyclic filter graphs and filters used in multi-mesh
161  // computation.
162 
163  if (this->solutions[i]->get_transform() == solutions_sub_idx[i])
164  this->solutions[i]->push_transform(son);
165  solutions_sub_idx[i] = this->solutions[i]->get_transform();
166  }
167  }
168 
169  template<typename Scalar>
171  {
173  for (int i = 0; i < this->solutions.size(); i++)
174  {
175  if (this->solutions[i]->get_transform() == solutions_sub_idx[i] && solutions_sub_idx[i] != 0)
176  this->solutions[i]->pop_transform();
177  solutions_sub_idx[i] = this->solutions[i]->get_transform();
178  }
179  }
180 
181  template<typename Scalar>
183  {
184  }
185 
186  template<typename Scalar>
187  SimpleFilter<Scalar>::SimpleFilter(std::vector<MeshFunctionSharedPtr<Scalar> > solutions, const std::vector<int> items) : Filter<Scalar>(solutions), items(items)
188  {
189  if (this->items.size() > 0)
190  Hermes::Helpers::check_length(solutions, items);
191  else
192  for (int i = 0; i < this->solutions.size(); i++)
193  this->items.push_back(H2D_FN_VAL);
194 
195  this->init();
196  this->init_components();
197  }
198 
199  template<typename Scalar>
200  SimpleFilter<Scalar>::~SimpleFilter()
201  {
202  }
203 
204  template<typename Scalar>
205  void SimpleFilter<Scalar>::init_components()
206  {
207  bool vec1 = false, vec2 = false;
208  for (int i = 0; i < this->solutions.size(); i++)
209  {
210  if (this->solutions[i]->get_num_components() > 1)
211  vec1 = true;
212 
213  if ((this->items[i] & H2D_FN_COMPONENT_0) && (this->items[i] & H2D_FN_COMPONENT_1))
214  vec2 = true;
215 
216  if (this->solutions[i]->get_num_components() == 1)
217  this->items[i] &= H2D_FN_COMPONENT_0;
218  }
219  this->num_components = (vec1 && vec2) ? 2 : 1;
220  }
221 
222  template<typename Scalar>
223  void SimpleFilter<Scalar>::precalculate(unsigned short order, unsigned short mask)
224  {
225 #ifdef H2D_USE_SECOND_DERIVATIVES
226  if (mask & (H2D_FN_DX | H2D_FN_DY | H2D_FN_DXX | H2D_FN_DYY | H2D_FN_DXY))
227 #else
228  if (mask & (H2D_FN_DX | H2D_FN_DY))
229 #endif
230  throw Hermes::Exceptions::Exception("SimpleFilter not defined for derivatives.");
231 
232  Quad2D* quad = this->quads[this->cur_quad];
233  unsigned char np = quad->get_num_points(order, this->element->get_mode());
234 
235  // precalculate all solutions
236  for (int i = 0; i < this->solutions.size(); i++)
237  this->solutions[i]->set_quad_order(order, this->items[i]);
238 
239  for (int j = 0; j < this->num_components; j++)
240  {
241  // obtain corresponding tables
242  std::vector<const Scalar*> tab;
243  for (int i = 0; i < this->solutions.size(); i++)
244  {
245  int a = 0, b = 0, mask = this->items[i];
246  if (mask >= 0x40) { a = 1; mask >>= 6; }
247  while (!(mask & 1)) { mask >>= 1; b++; }
248  tab.push_back(const_cast<Scalar*>(this->solutions[i]->get_values(this->num_components == 1 ? a : j, b)));
249  if (tab[i] == nullptr)
250  throw Hermes::Exceptions::Exception("Value of 'item%d' is incorrect in filter definition.", i + 1);
251  }
252 
253  // apply the filter
254  filter_fn(np, tab, this->values[j][0]);
255  }
256  this->values_valid = true;
257  }
258 
259  template<typename Scalar>
260  Func<Scalar>* SimpleFilter<Scalar>::get_pt_value(double x, double y, bool use_MeshHashGrid, Element* e)
261  {
262  std::vector<Scalar> val;
263  for (int i = 0; i < this->solutions.size(); i++)
264  val.push_back(this->solutions[i]->get_pt_value(x, y, use_MeshHashGrid, e)->val[0]);
265 
266  Func<Scalar>* toReturn = new Func<Scalar>(1, 1);
267 
268  Scalar result;
269  std::vector<const Scalar*> values;
270  for (int i = 0; i < this->solutions.size(); i++)
271  values.push_back(&val[i]);
272 
273  // apply the filter
274  filter_fn(1, values, &result);
275 
276  toReturn->val[0] = result;
277  return toReturn;
278  }
279 
280  ComplexFilter::ComplexFilter() : Filter<double>()
281  {
282  this->unimesh = false;
283  }
284 
285  ComplexFilter::ComplexFilter(MeshFunctionSharedPtr<std::complex<double> > solution, int item) : Filter<double>()
286  {
287  this->unimesh = false;
288  this->sln_complex = solution;
289  this->num_components = solution->get_num_components();
290  this->mesh = solution->get_mesh();
291  this->solutions_sub_idx.push_back(0);
292  set_quad_2d(&g_quad_2d_std);
293  }
294 
295  ComplexFilter::~ComplexFilter()
296  {
297  this->free();
298  }
299 
301  {
302  }
303 
305  {
307  this->sln_complex->set_quad_2d(quad_2d);
308  }
309 
311  {
313  this->sln_complex->set_active_element(e);
314  this->solutions_sub_idx[0] = 0;
315 
316  // fixme
317  this->order = 20;
318  }
319 
321  {
323  this->sln_complex->push_transform(son);
324  }
325 
327  {
329  this->sln_complex->pop_transform();
330  }
331 
332  void ComplexFilter::precalculate(unsigned short order, unsigned short mask)
333  {
334 #ifdef H2D_USE_SECOND_DERIVATIVES
335  if (mask & (H2D_FN_DX | H2D_FN_DY | H2D_FN_DXX | H2D_FN_DYY | H2D_FN_DXY))
336 #else
337  if (mask & (H2D_FN_DX | H2D_FN_DY))
338 #endif
339  throw Hermes::Exceptions::Exception("Filter not defined for derivatives.");
340 
341  Quad2D* quad = this->quads[this->cur_quad];
342  unsigned char np = quad->get_num_points(order, this->element->get_mode());
343 
344  this->sln_complex->set_quad_order(order, H2D_FN_VAL);
345 
346  // obtain corresponding tables
347  filter_fn(np, const_cast<std::complex<double>*>(this->sln_complex->get_values(0, 0)), this->values[0][0]);
348  if (num_components > 1)
349  filter_fn(np, const_cast<std::complex<double>*>(this->sln_complex->get_values(1, 0)), this->values[1][0]);
350  this->values_valid = true;
351  }
352 
353  Func<double>* ComplexFilter::get_pt_value(double x, double y, bool use_MeshHashGrid, Element* e)
354  {
355  Func<std::complex<double> >* val = this->sln_complex->get_pt_value(x, y, use_MeshHashGrid, e);
356 
357  Func<double>* toReturn = new Func<double>(1, this->num_components);
358 
359  double result;
360 
361  // apply the filter
362  filter_fn(1, val->val, &result);
363 
364  if (this->num_components == 1)
365  {
366  toReturn->val[0] = result;
367  toReturn->dx[0] = 0.0;
368  toReturn->dy[0] = 0.0;
369  }
370  else
371  {
372  this->warn("ComplexFilter only implemented for scalar functions.");
373  }
374  return toReturn;
375  }
376 
377  template<typename Scalar>
379  {
380  }
381 
382  template<typename Scalar>
383  DXDYFilter<Scalar>::DXDYFilter(std::vector<MeshFunctionSharedPtr<Scalar> > solutions) : Filter<Scalar>(solutions)
384  {
385  init_components();
386  }
387 
388  template<typename Scalar>
389  DXDYFilter<Scalar>::~DXDYFilter()
390  {
391  }
392 
393  template<typename Scalar>
394  void DXDYFilter<Scalar>::init_components()
395  {
396  this->num_components = this->solutions[0]->get_num_components();
397  for (int i = 1; i < this->solutions.size(); i++)
398  if (this->solutions[i]->get_num_components() != this->num_components)
399  throw Hermes::Exceptions::Exception("Filter: Solutions do not have the same number of components!");
400  }
401 
402  template<typename Scalar>
403  void DXDYFilter<Scalar>::init(std::vector<MeshFunctionSharedPtr<Scalar> > solutions)
404  {
406  init_components();
407  }
408 
409  template<typename Scalar>
410  void DXDYFilter<Scalar>::precalculate(unsigned short order, unsigned short mask)
411  {
412  Quad2D* quad = this->quads[this->cur_quad];
413  unsigned char np = quad->get_num_points(order, this->element->get_mode());
414 
415  // precalculate all solutions
416  for (int i = 0; i < this->solutions.size(); i++)
417  this->solutions[i]->set_quad_order(order, H2D_FN_DEFAULT);
418 
419  for (int j = 0; j < this->num_components; j++)
420  {
421  // obtain solution tables
422  double *x, *y;
423  const Scalar *val[H2D_MAX_COMPONENTS], *dx[H2D_MAX_COMPONENTS], *dy[H2D_MAX_COMPONENTS];
424  x = this->solutions[0]->get_refmap()->get_phys_x(order);
425  y = this->solutions[0]->get_refmap()->get_phys_y(order);
426 
427  for (int i = 0; i < this->solutions.size(); i++)
428  {
429  val[i] = this->solutions[i]->get_fn_values(j);
430  dx[i] = this->solutions[i]->get_dx_values(j);
431  dy[i] = this->solutions[i]->get_dy_values(j);
432  }
433 
434  std::vector<const Scalar *> values_vector;
435  std::vector<const Scalar *> dx_vector;
436  std::vector<const Scalar *> dy_vector;
437 
438  for (int i = 0; i < this->solutions.size(); i++)
439  {
440  values_vector.push_back(val[i]);
441  dx_vector.push_back(dx[i]);
442  dy_vector.push_back(dy[i]);
443  }
444 
445  // apply the filter
446  filter_fn(np, x, y, values_vector, dx_vector, dy_vector, this->values[j][0], this->values[j][1], this->values[j][2]);
447  }
448 
449  this->values_valid = true;
450  }
451 
452  template<typename Scalar>
453  Func<Scalar>* DXDYFilter<Scalar>::get_pt_value(double x, double y, bool use_MeshHashGrid, Element* e)
454  {
455  this->warn("DXDYFilter<Scalar>::get_pt_value not implemented.");
456  return 0;
457  }
458 
459  template<typename Scalar>
460  void MagFilter<Scalar>::filter_fn(int n, const std::vector<const Scalar*>& values, Scalar* result)
461  {
462  for (int i = 0; i < n; i++)
463 
464  {
465  result[i] = 0;
466  for (unsigned int j = 0; j < values.size(); j++)
467  result[i] += sqr(values[j][i]);
468  result[i] = sqrt(result[i]);
469  }
470  };
471 
472  template<typename Scalar>
473  MagFilter<Scalar>::MagFilter(std::vector<MeshFunctionSharedPtr<Scalar> > solutions, std::vector<int> items) : SimpleFilter<Scalar>(solutions, items)
474  {
476  };
477 
478  template<typename Scalar>
480  : SimpleFilter<Scalar>()
481  {
482  this->solutions = { sln1, sln1 };
483  this->items = { item1 & H2D_FN_COMPONENT_0, item1 & H2D_FN_COMPONENT_1 };
484  if (sln1->get_num_components() < 2)
485  throw Hermes::Exceptions::Exception("The single-argument constructor is intended for vector-valued solutions.");
487  };
488 
489  template<typename Scalar>
491  {
492  };
493 
494  template<typename Scalar>
496  {
497  std::vector<MeshFunctionSharedPtr<Scalar> > slns;
498  std::vector<int> items;
499  for (int i = 0; i < this->solutions.size(); i++)
500  {
501  slns.push_back(this->solutions[i]->clone());
502  items.push_back(this->items[i]);
503  }
504  MagFilter<Scalar>* filter = new MagFilter<Scalar>(slns, items);
505  return filter;
506  }
507 
508  void TopValFilter::filter_fn(int n, const std::vector<const double*>& values, double* result)
509  {
510  for (int i = 0; i < n; i++)
511  {
512  result[i] = 0;
513  for (unsigned int j = 0; j < values.size(); j++)
514  if (values[j][i] > limits[j])
515  result[i] = limits[j];
516  else
517  result[i] = values[j][i];
518  }
519  };
520 
521  TopValFilter::TopValFilter(std::vector<MeshFunctionSharedPtr<double> > solutions, std::vector<double> limits, std::vector<int> items) : SimpleFilter<double>(solutions, items), limits(limits)
522  {
523  };
524 
525  TopValFilter::TopValFilter(MeshFunctionSharedPtr<double> sln, double limit, int item)
526  : SimpleFilter<double>()
527  {
528  this->limits.push_back(limit);
529  this->solutions.push_back(sln);
530  this->items.push_back(item);
532  };
533 
534  TopValFilter::~TopValFilter()
535  {
536  }
537 
539  {
540  std::vector<MeshFunctionSharedPtr<double> > slns;
541  std::vector<int> items;
542  for (int i = 0; i < this->solutions.size(); i++)
543  {
544  slns.push_back(this->solutions[i]->clone());
545  items.push_back(this->items[i]);
546  }
547  TopValFilter* filter = new TopValFilter(slns, limits, items);
548  return filter;
549  }
550 
551  void BottomValFilter::filter_fn(int n, const std::vector<const double*>& values, double* result)
552  {
553  for (int i = 0; i < n; i++)
554  {
555  result[i] = 0;
556  for (unsigned int j = 0; j < values.size(); j++)
557  if (values[j][i] < limits[j])
558  result[i] = limits[j];
559  else
560  result[i] = values[j][i];
561  }
562  };
563 
564  BottomValFilter::BottomValFilter(std::vector<MeshFunctionSharedPtr<double> > solutions, std::vector<double> limits, std::vector<int> items) : SimpleFilter<double>(solutions, items), limits(limits)
565  {
566  };
567 
568  BottomValFilter::BottomValFilter(MeshFunctionSharedPtr<double> sln, double limit, int item)
569  : SimpleFilter<double>()
570  {
571  this->limits.push_back(limit);
572  this->solutions.push_back(sln);
573  this->items.push_back(item);
575  };
576 
577  BottomValFilter::~BottomValFilter()
578  {
579  }
580 
582  {
583  std::vector<MeshFunctionSharedPtr<double> > slns;
584  std::vector<int> items;
585  for (int i = 0; i < this->solutions.size(); i++)
586  {
587  slns.push_back(this->solutions[i]->clone());
588  items.push_back(this->items[i]);
589  }
590  BottomValFilter* filter = new BottomValFilter(slns, limits, items);
591  return filter;
592  }
593 
594  void ValFilter::filter_fn(int n, const std::vector<const double*>& values, double* result)
595  {
596  for (int i = 0; i < n; i++)
597  {
598  result[i] = 0;
599  for (unsigned int j = 0; j < values.size(); j++)
600  if (values[j][i] < low_limits[j])
601  result[i] = low_limits[j];
602  else
603  if (values[j][i] > high_limits[j])
604  result[i] = high_limits[j];
605  else
606  result[i] = values[j][i];
607  }
608  };
609 
610  ValFilter::ValFilter(std::vector<MeshFunctionSharedPtr<double> > solutions, std::vector<double> low_limits, std::vector<double> high_limits, std::vector<int> items) : SimpleFilter<double>(solutions, items), low_limits(low_limits), high_limits(high_limits)
611  {
612  };
613 
614  ValFilter::ValFilter(MeshFunctionSharedPtr<double> sln, double low_limit, double high_limit, int item)
615  : SimpleFilter<double>()
616  {
617  this->low_limits.push_back(low_limit);
618  this->high_limits.push_back(high_limit);
619  this->solutions.push_back(sln);
620  this->items.push_back(item);
622  };
623 
624  ValFilter::~ValFilter()
625  {
626  }
627 
629  {
630  std::vector<MeshFunctionSharedPtr<double> > slns;
631  std::vector<int> items;
632  for (int i = 0; i < this->solutions.size(); i++)
633  {
634  slns.push_back(this->solutions[i]->clone());
635  items.push_back(this->items[i]);
636  }
637  ValFilter* filter = new ValFilter(slns, low_limits, high_limits, items);
638  return filter;
639  }
640 
641  template<typename Scalar>
642  void DiffFilter<Scalar>::filter_fn(int n, const std::vector<const Scalar*>& values, Scalar* result)
643  {
644  for (int i = 0; i < n; i++) result[i] = values[0][i] - values.at(1)[i];
645  };
646 
647  template<typename Scalar>
648  DiffFilter<Scalar>::DiffFilter(std::vector<MeshFunctionSharedPtr<Scalar> > solutions, std::vector<int> items) : SimpleFilter<Scalar>(solutions, items) {}
649 
650  template<typename Scalar>
651  DiffFilter<Scalar>::~DiffFilter()
652  {
653  }
654 
655  template<typename Scalar>
657  {
658  std::vector<MeshFunctionSharedPtr<Scalar> > slns;
659  std::vector<int> items;
660  for (int i = 0; i < this->solutions.size(); i++)
661  {
662  slns.push_back(this->solutions[i]->clone());
663  items.push_back(this->items[i]);
664  }
665  DiffFilter* filter = new DiffFilter<Scalar>(slns, items);
666  return filter;
667  }
668 
669  template<typename Scalar>
670  void SumFilter<Scalar>::filter_fn(int n, const std::vector<const Scalar*>& values, Scalar* result)
671  {
672  for (int i = 0; i < n; i++)
673  {
674  result[i] = 0;
675  for (unsigned int j = 0; j < values.size(); j++)
676  result[i] += values[j][i];
677  }
678  };
679 
680  template<typename Scalar>
681  SumFilter<Scalar>::SumFilter(std::vector<MeshFunctionSharedPtr<Scalar> > solutions, std::vector<int> items) : SimpleFilter<Scalar>(solutions, items) {}
682 
683  template<typename Scalar>
684  SumFilter<Scalar>::~SumFilter()
685  {
686  }
687 
688  template<typename Scalar>
690  {
691  std::vector<MeshFunctionSharedPtr<Scalar> > slns;
692  std::vector<int> items;
693  for (int i = 0; i < this->solutions.size(); i++)
694  {
695  slns.push_back(this->solutions[i]->clone());
696  items.push_back(this->items[i]);
697  }
698  SumFilter<Scalar>* filter = new SumFilter<Scalar>(slns, items);
699  return filter;
700  }
701 
702  template<>
703  void SquareFilter<double>::filter_fn(int n, const std::vector<const double *>& v1, double* result)
704  {
705  for (int i = 0; i < n; i++)
706  result[i] = sqr(v1[0][i]);
707  };
708 
709  template<>
710  void SquareFilter<std::complex<double> >::filter_fn(int n, const std::vector<const std::complex<double> *>& v1, std::complex<double> * result)
711  {
712  for (int i = 0; i < n; i++)
713  result[i] = std::norm(v1[0][i]);
714  };
715 
716  template<typename Scalar>
717  SquareFilter<Scalar>::SquareFilter(std::vector<MeshFunctionSharedPtr<Scalar> > solutions, std::vector<int> items)
718  : SimpleFilter<Scalar>(solutions, items)
719  {
720  if (solutions.size() > 1)
721  throw Hermes::Exceptions::Exception("SquareFilter only supports one MeshFunction.");
722  };
723 
724  template<typename Scalar>
725  SquareFilter<Scalar>::~SquareFilter()
726  {
727  }
728 
729  template<typename Scalar>
731  {
732  std::vector<MeshFunctionSharedPtr<Scalar> > slns;
733  std::vector<int> items;
734  for (int i = 0; i < this->solutions.size(); i++)
735  {
736  slns.push_back(this->solutions[i]->clone());
737  items.push_back(this->items[i]);
738  }
739  SquareFilter<Scalar>* filter = new SquareFilter<Scalar>(slns, items);
740  return filter;
741  }
742 
743  void AbsFilter::filter_fn(int n, const std::vector<const double*>& v1, double * result)
744  {
745  for (int i = 0; i < n; i++)
746  result[i] = std::abs(v1[0][i]);
747  };
748 
749  AbsFilter::AbsFilter(std::vector<MeshFunctionSharedPtr<double> > solutions, std::vector<int> items)
750  : SimpleFilter<double>(solutions, items)
751  {
752  if (solutions.size() > 1)
753  throw Hermes::Exceptions::Exception("AbsFilter only supports one MeshFunction.");
754  };
755 
756  AbsFilter::AbsFilter(MeshFunctionSharedPtr<double> solution)
757  : SimpleFilter<double>()
758  {
759  this->solutions.push_back(solution);
760 
761  this->items.push_back(H2D_FN_VAL);
762 
763  this->init();
764  init_components();
765  };
766 
767  AbsFilter::~AbsFilter()
768  {
769  }
770 
772  {
773  std::vector<MeshFunctionSharedPtr<double> > slns;
774  std::vector<int> items;
775  for (int i = 0; i < this->solutions.size(); i++)
776  {
777  slns.push_back(this->solutions[i]->clone());
778  items.push_back(this->items[i]);
779  }
780  AbsFilter* filter = new AbsFilter(slns, items);
781  return filter;
782  }
783 
784  void RealFilter::filter_fn(int n, const std::complex<double>* values, double* result)
785  {
786  for (int i = 0; i < n; i++)
787  result[i] = values[i].real();
788  };
789 
791  {
792  RealFilter* filter = new RealFilter(this->sln_complex->clone(), this->item);
793  return filter;
794  }
795 
796  RealFilter::RealFilter()
797  : ComplexFilter()
798  {
799  };
800 
801  RealFilter::RealFilter(MeshFunctionSharedPtr<std::complex<double> > solution, int item)
802  : ComplexFilter(solution, item)
803  {
804  };
805 
806  RealFilter::~RealFilter()
807  {
808  }
809 
810  void ImagFilter::filter_fn(int n, const std::complex<double>* values, double* result)
811  {
812  for (int i = 0; i < n; i++)
813  result[i] = values[i].imag();
814  };
815 
816  ImagFilter::ImagFilter(MeshFunctionSharedPtr<std::complex<double> > solution, int item)
817  : ComplexFilter(solution, item)
818  {
819  };
820 
821  ImagFilter::~ImagFilter()
822  {
823  }
824 
826  {
827  ImagFilter* filter = new ImagFilter(this->sln_complex->clone(), this->item);
828  return filter;
829  }
830 
831  void ComplexAbsFilter::filter_fn(int n, const std::complex<double>* values, double* result)
832  {
833  for (int i = 0; i < n; i++)
834  result[i] = sqrt(sqr(values[i].real()) + sqr(values[i].imag()));
835  };
836 
838  {
839  ComplexAbsFilter* filter = new ComplexAbsFilter(this->sln_complex->clone(), this->item);
840  return filter;
841  }
842 
843  ComplexAbsFilter::ComplexAbsFilter(MeshFunctionSharedPtr<std::complex<double> > solution, int item)
844  : ComplexFilter(solution, item)
845  {
846  };
847 
848  ComplexAbsFilter::~ComplexAbsFilter()
849  {
850  }
851 
852  void AngleFilter::filter_fn(int n, const std::vector<const std::complex<double>*>& v1, double* result)
853  {
854  for (int i = 0; i < n; i++)
855  result[i] = atan2(v1[0][i].imag(), v1[0][i].real());
856  };
857 
858  AngleFilter::AngleFilter(std::vector<MeshFunctionSharedPtr<std::complex<double> > > solutions, std::vector<int> items)
859  : SimpleFilter<std::complex<double> >(solutions, items)
860  {
861  if (solutions.size() > 1)
862  throw Hermes::Exceptions::Exception("RealFilter only supports one MeshFunction.");
863  };
864 
865  AngleFilter::~AngleFilter()
866  {
867  }
868 
869  void VonMisesFilter::precalculate(unsigned short order, unsigned short mask)
870  {
871 #ifdef H2D_USE_SECOND_DERIVATIVES
872  if (mask & (H2D_FN_DX | H2D_FN_DY | H2D_FN_DXX | H2D_FN_DYY | H2D_FN_DXY))
873 #else
874  if (mask & (H2D_FN_DX | H2D_FN_DY))
875 #endif
876  throw Hermes::Exceptions::Exception("VonMisesFilter not defined for derivatives.");
877 
878  Quad2D* quad = this->quads[this->cur_quad];
879  unsigned char np = quad->get_num_points(order, this->element->get_mode());
880 
881  this->solutions[0]->set_quad_order(order, H2D_FN_VAL | H2D_FN_DX | H2D_FN_DY);
882  this->solutions[1]->set_quad_order(order, H2D_FN_DX | H2D_FN_DY);
883 
884  const double *dudx = this->solutions[0]->get_dx_values();
885  const double *dudy = this->solutions[0]->get_dy_values();
886  const double *dvdx = this->solutions[1]->get_dx_values();
887  const double *dvdy = this->solutions[1]->get_dy_values();
888  const double *uval = this->solutions[0]->get_fn_values();
889  update_refmap();
890  double *x = refmap.get_phys_x(order);
891 
892  for (int i = 0; i < np; i++)
893  {
894  // stress tensor
895  double tz = lambda*(dudx[i] + dvdy[i]);
896  double tx = tz + 2 * mu*dudx[i];
897  double ty = tz + 2 * mu*dvdy[i];
898  if (cyl) tz += 2 * mu*uval[i] / x[i];
899  double txy = mu*(dudy[i] + dvdx[i]);
900 
901  // Von Mises stress
902  this->values[0][0][i] = 1.0 / sqrt(2.0) * sqrt(sqr(tx - ty) + sqr(ty - tz) + sqr(tz - tx) + 6 * sqr(txy));
903  }
904  this->values_valid = true;
905  }
906 
907  Func<double>* VonMisesFilter::get_pt_value(double x, double y, bool use_MeshHashGrid, Element* e)
908  {
909  this->warn("VonMisesFilter<Scalar>::get_pt_value not implemented.");
910  return 0;
911  }
912 
913  VonMisesFilter::VonMisesFilter(std::vector<MeshFunctionSharedPtr<double> > solutions, double lambda, double mu,
914  int cyl, int item1, int item2)
915  : Filter<double>(solutions)
916  {
917  this->mu = mu;
918  this->lambda = lambda;
919  this->cyl = cyl;
920  this->item1 = item1;
921  this->item2 = item2;
922  }
923 
924  VonMisesFilter::~VonMisesFilter()
925  {
926  }
927 
929  {
930  std::vector<MeshFunctionSharedPtr<double> > slns;
931  for (int i = 0; i < this->solutions.size(); i++)
932  slns.push_back(this->solutions[i]->clone());
933  VonMisesFilter* filter = new VonMisesFilter(slns, lambda, mu, cyl, item1, item2);
934  return filter;
935  }
936 
937  template<typename Scalar>
938  void LinearFilter<Scalar>::precalculate(unsigned short order, unsigned short mask)
939  {
940  Quad2D* quad = this->quads[this->cur_quad];
941  unsigned char np = quad->get_num_points(order, this->element->get_mode());
942  struct Filter<Scalar>::Node* node = this->new_node(H2D_FN_DEFAULT, np);
943 
944  // precalculate all solutions
945  for (int i = 0; i < this->solutions.size(); i++)
946  this->solutions[i]->set_quad_order(order);
947 
948  for (int j = 0; j < this->num_components; j++)
949  {
950  // obtain solution tables
951  std::vector<Scalar*> val, dx, dy;
952  for (int i = 0; i < this->solutions.size(); i++)
953  {
954  val.push_back(this->solutions[i]->get_fn_values(j));
955  dx.push_back(this->solutions[i]->get_dx_values(j));
956  dy.push_back(this->solutions[i]->get_dy_values(j));
957  }
958  if (this->solutions.size() == 2)
959  {
960  for (int i = 0; i < np; i++)
961  {
962  node->values[j][0][i] = tau_frac * (val[1][i] - val[0][i]) + val[1][i];
963  node->values[j][1][i] = tau_frac * (dx[1][i] - dx[0][i]) + dx[1][i];
964  node->values[j][2][i] = tau_frac * (dy[1][i] - dy[0][i]) + dy[1][i];
965  }
966  }
967  else
968  {
969  for (int i = 0; i < np; i++)
970  {
971  node->values[j][0][i] = val[0][i];
972  node->values[j][1][i] = dx[0][i];
973  node->values[j][2][i] = dy[0][i];
974  }
975  }
976  }
977  this->values_valid = true;
978  }
979 
980  template<typename Scalar>
981  Func<Scalar>* LinearFilter<Scalar>::get_pt_value(double x, double y, bool use_MeshHashGrid, Element* e)
982  {
983  this->warn("LinearFilter<Scalar>::get_pt_value not implemented.");
984  return 0;
985  }
986 
987  template<typename Scalar>
989  {
990  init_components();
991  }
992 
993  template<typename Scalar>
994  LinearFilter<Scalar>::LinearFilter(MeshFunctionSharedPtr<Scalar> older, MeshFunctionSharedPtr<Scalar> old, double tau_frac)
995  : Filter<Scalar>({ older, old })
996  {
997  this->tau_frac = tau_frac;
998  init_components();
999  }
1000 
1001  template<typename Scalar>
1002  LinearFilter<Scalar>::~LinearFilter()
1003  {
1004  }
1005 
1006  template<typename Scalar>
1007  void LinearFilter<Scalar>::init_components()
1008  {
1009  this->num_components = this->solutions[0]->get_num_components();
1010  for (int i = 1; i < this->solutions.size(); i++)
1011  if (this->solutions[i]->get_num_components() != this->num_components)
1012  throw Hermes::Exceptions::Exception("Filter: Solutions do not have the same number of components!");
1013  }
1014 
1015  template<typename Scalar>
1017  {
1019 
1020  this->order = 0;
1021  for (int i = 0; i < this->solutions.size(); i++)
1022  {
1023  int o = this->solutions[i]->get_fn_order();
1024  if (o > this->order) this->order = o;
1025  }
1026  }
1027 
1028  template class HERMES_API Filter < double > ;
1029  template class HERMES_API Filter < std::complex<double> > ;
1030  template class HERMES_API SimpleFilter < double > ;
1031  template class HERMES_API SimpleFilter < std::complex<double> > ;
1032  template class HERMES_API DXDYFilter < double > ;
1033  template class HERMES_API DXDYFilter < std::complex<double> > ;
1034  template class HERMES_API MagFilter < double > ;
1035  template class HERMES_API MagFilter < std::complex<double> > ;
1036  template class HERMES_API DiffFilter < double > ;
1037  template class HERMES_API DiffFilter < std::complex<double> > ;
1038  template class HERMES_API SumFilter < double > ;
1039  template class HERMES_API SumFilter < std::complex<double> > ;
1040  template class HERMES_API SquareFilter < double > ;
1041  template class HERMES_API SquareFilter < std::complex<double> > ;
1042  }
1043 }
virtual void init()
Internal.
Definition: filter.cpp:38
virtual MeshFunction< double > * clone() const
Definition: filter.cpp:581
Definition: adapt.h:24
virtual void push_transform(int son)
Definition: function.cpp:115
int id
element id number
Definition: element.h:112
Calculates the Von Mises stress.
Definition: filter.h:343
Stores one element of a mesh.
Definition: element.h:107
Represents a finite element mesh. Typical usage: MeshSharedPtr mesh; Hermes::Hermes2D::MeshReaderH2DX...
Definition: mesh.h:61
virtual void set_active_element(Element *e)
virtual void set_active_element(Element *e)
Definition: filter.cpp:1016
TopValFilter takes functions and puts a threshold on their highest values.
Definition: filter.h:188
virtual void set_active_element(Element *e)
Definition: filter.cpp:104
virtual MeshFunction< double > * clone() const
Definition: filter.cpp:837
virtual void set_quad_2d(Quad2D *quad_2d)
Represents a function defined on a mesh.
Definition: mesh_function.h:56
virtual void set_quad_2d(Quad2D *quad_2d)
Definition: filter.cpp:95
virtual void push_transform(int son)
Definition: filter.cpp:149
virtual MeshFunction< double > * clone() const
Definition: filter.cpp:628
virtual MeshFunction< Scalar > * clone() const
Definition: filter.cpp:689
virtual void precalculate(unsigned short order, unsigned short mask)
precalculates the current function at the current integration points.
Definition: filter.cpp:223
virtual void reinit()
Internal.
Definition: filter.cpp:142
virtual MeshFunction< double > * clone() const
Definition: filter.cpp:771
double * get_phys_x(int order)
Definition: refmap.cpp:38
virtual MeshFunction< Scalar > * clone() const
Definition: filter.cpp:656
virtual void set_quad_2d(Quad2D *quad_2d)
Definition: filter.cpp:304
Stores one node of a mesh.
Definition: element.h:45
virtual void pop_transform()
Definition: filter.cpp:170
virtual const Scalar * get_dx_values(int component=0) const
Returns the x partial derivative.
Definition: function.cpp:175
virtual const Scalar * get_dy_values(int component=0) const
Returns the y partial derivative.
Definition: function.cpp:182
virtual void precalculate(unsigned short order, unsigned short mask)
precalculates the current function at the current integration points.
Definition: filter.cpp:938
BottomValFilter takes functions and puts a threshold on their lowest values.
Definition: filter.h:204
Calculates the magnitude of a vector function.
Definition: filter.h:173
virtual MeshFunction< Scalar > * clone() const
Definition: filter.cpp:730
virtual void push_transform(int son)
Definition: filter.cpp:320
Calculates absolute value of a real solution.
Definition: filter.h:276
virtual void set_active_element(Element *e)
Definition: filter.cpp:310
ValFilter takes functions and puts a threshold on their lowest AND highest values.
Definition: filter.h:220
const int H2D_FN_DEFAULT
default precalculation mask
Definition: function.h:57
virtual void pop_transform()
Definition: function.cpp:122
Calculates the square of a function.
Definition: filter.h:264
virtual MeshFunction< double > * clone() const
Definition: filter.cpp:928
double values[H2D_MAX_SOLUTION_COMPONENTS][H2D_NUM_FUNCTION_VALUES][H2D_MAX_INTEGRATION_POINTS_COUNT]
The data.
Definition: function.h:203
virtual MeshFunction< double > * clone() const
Definition: filter.cpp:538
virtual void free()
Frees all precalculated tables.
Definition: filter.cpp:131
Calculates the difference of two functions.
Definition: filter.h:238
virtual void precalculate(unsigned short order, unsigned short mask)
precalculates the current function at the current integration points.
Definition: filter.cpp:869
Computes the absolute value of a complex solution.
Definition: filter.h:316
virtual Func< Scalar > * get_pt_value(double x, double y, bool use_MeshHashGrid=false, Element *e=nullptr)
Return the value at the coordinates x,y.
Definition: filter.cpp:981
Calculates the sum of two functions.
Definition: filter.h:251
virtual Func< double > * get_pt_value(double x, double y, bool use_MeshHashGrid=false, Element *e=nullptr)
Return the value at the coordinates x,y.
Definition: filter.cpp:353
Element * element
The active element.
bool values_valid
Flag that the data are not 'dirty'.
Definition: function.h:205
Calculated function values (from the class Function) on an element for assembling.
Definition: forms.h:214
VonMisesFilter(std::vector< MeshFunctionSharedPtr< double > > solutions, double lambda, double mu, int cyl=0, int item1=H2D_FN_VAL, int item2=H2D_FN_VAL)
Definition: filter.cpp:913
const int H2D_FN_VAL
Both components are usually requested together...
Definition: function.h:53
unsigned char num_components
Number of vector components.
Definition: function.h:218
Removes the imaginary part from a function.
Definition: filter.h:289
virtual Func< Scalar > * get_pt_value(double x, double y, bool use_MeshHashGrid=false, Element *e=nullptr)
Return the value at the coordinates x,y.
Definition: filter.cpp:260
virtual const Scalar * get_fn_values(int component=0) const
Returns function values.
Definition: function.cpp:168
virtual MeshFunction< Scalar > * clone() const
Definition: filter.cpp:495
virtual void free()
Frees all precalculated tables.
Definition: filter.cpp:300
virtual MeshFunction< double > * clone() const
Definition: filter.cpp:825
Quad2D * quads[H2D_MAX_QUADRATURES]
List of available quadratures.
Definition: function.h:226
int cur_quad
Active quadrature (index into 'quads')
Definition: function.h:228
virtual void precalculate(unsigned short order, unsigned short mask)
precalculates the current function at the current integration points.
Definition: filter.cpp:410
Represents the solution of a PDE.
Definition: api2d.h:35
virtual void precalculate(unsigned short order, unsigned short mask)
precalculates the current function at the current integration points.
Definition: filter.cpp:332
virtual MeshFunction< double > * clone() const
Definition: filter.cpp:790
virtual void pop_transform()
Definition: filter.cpp:326
int order
Current function polynomial order.
Definition: function.h:215
virtual Func< double > * get_pt_value(double x, double y, bool use_MeshHashGrid=false, Element *e=nullptr)
Return the value at the coordinates x,y.
Definition: filter.cpp:907
virtual Func< Scalar > * get_pt_value(double x, double y, bool use_MeshHashGrid=false, Element *e=nullptr)
Return the value at the coordinates x,y.
Definition: filter.cpp:453