Hermes2D  2.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(MeshFunction<Scalar>** solutions, int num) : MeshFunction<Scalar>()
33  {
34  this->num = num;
35  if(num > 10)
36  throw Hermes::Exceptions::Exception("Attempt to create an instance of Filter with more than 10 MeshFunctions.");
37  for(int i = 0; i < this->num; i++)
38  this->sln[i] = solutions[i];
39  this->init();
40  }
41 
42  template<typename Scalar>
43  Filter<Scalar>::Filter(const Hermes::vector<MeshFunction<Scalar>*>& solutions) : MeshFunction<Scalar>()
44  {
45  this->num = solutions.size();
46  if(num > 10)
47  throw Hermes::Exceptions::Exception("Attempt to create an instance of Filter with more than 10 MeshFunctions.");
48  for(int i = 0; i < this->num; i++)
49  this->sln[i] = solutions.at(i);
50  this->init();
51  }
52 
53  template<typename Scalar>
54  Filter<Scalar>::Filter(const Hermes::vector<Solution<Scalar>*>& solutions) : MeshFunction<Scalar>()
55  {
56  this->num = solutions.size();
57  if(num > 10)
58  throw Hermes::Exceptions::Exception("Attempt to create an instance of Filter with more than 10 MeshFunctions.");
59  for(int i = 0; i < this->num; i++)
60  this->sln[i] = solutions.at(i);
61  this->init();
62  }
63 
64  template<typename Scalar>
65  void Filter<Scalar>::init(const Hermes::vector<MeshFunction<Scalar>*>& solutions)
66  {
67  this->num = solutions.size();
68  if(num > 10)
69  throw Hermes::Exceptions::Exception("Attempt to create an instance of Filter with more than 10 MeshFunctions.");
70  for(int i = 0; i < this->num; i++)
71  this->sln[i] = solutions.at(i);
72  this->init();
73  }
74 
75  template<typename Scalar>
76  void Filter<Scalar>::setDeleteSolutions()
77  {
78  this->deleteSolutions = true;
79  }
80 
81  template<typename Scalar>
82  void Filter<Scalar>::init()
83  {
84  // construct the union mesh, if necessary
85  const Mesh* meshes[H2D_MAX_COMPONENTS];
86  for(int i = 0; i < this->num; i++)
87  meshes[i] = this->sln[i]->get_mesh();
88  this->mesh = meshes[0];
89  unimesh = false;
90 
91  for (int i = 1; i < num; i++)
92 
93  {
94  if(meshes[i] == NULL)
95  {
96  this->warn("You may be initializing a Filter with Solution that is missing a Mesh.");
97  throw Hermes::Exceptions::Exception("this->meshes[%d] is NULL in Filter<Scalar>::init().", i);
98  }
99  if(meshes[i]->get_seq() != this->mesh->get_seq())
100  {
101  unimesh = true;
102  break;
103  }
104  }
105 
106  if(unimesh)
107  {
108  Traverse trav;
109  trav.begin(num, meshes);
110  this->mesh = new Mesh;
111  unidata = trav.construct_union_mesh(const_cast<Hermes2D::Mesh*>(this->mesh));
112  trav.finish();
113  }
114 
115  // misc init
116  this->num_components = 1;
117  this->order = 0;
118 
119  memset(sln_sub, 0, sizeof(sln_sub));
120  set_quad_2d(&g_quad_2d_std);
121 
122  this->deleteSolutions = false;
123  }
124 
125  template<typename Scalar>
126  Filter<Scalar>::~Filter()
127  {
128  this->free();
129  if(this->deleteSolutions)
130  {
131  for(int i = 0; i < this->num; i++)
132  delete this->sln[i];
133  }
134  }
135 
136  template<typename Scalar>
137  void Filter<Scalar>::set_quad_2d(Quad2D* quad_2d)
138  {
139  MeshFunction<Scalar>::set_quad_2d(quad_2d);
140  for (int i = 0; i < num; i++)
141  this->sln[i]->set_quad_2d(quad_2d); // nodup
142  }
143 
144  template<typename Scalar>
145  void Filter<Scalar>::set_active_element(Element* e)
146  {
147  MeshFunction<Scalar>::set_active_element(e);
148  if(!unimesh)
149  {
150  for (int i = 0; i < num; i++)
151  this->sln[i]->set_active_element(e); // nodup
152  memset(sln_sub, 0, sizeof(sln_sub));
153  }
154  else
155  {
156  for (int i = 0; i < num; i++)
157  {
158  this->sln[i]->set_active_element(unidata[i][e->id].e);
159  this->sln[i]->set_transform(unidata[i][e->id].idx);
160  sln_sub[i] = this->sln[i]->get_transform();
161  }
162  }
163 
164  for(typename std::map<uint64_t, LightArray<struct Filter<Scalar>::Node*>*>::iterator it = tables[this->cur_quad].begin(); it != tables[this->cur_quad].end(); it++)
165  {
166  for(unsigned int l = 0; l < it->second->get_size(); l++)
167  if(it->second->present(l))
168  ::free(it->second->get(l));
169  delete it->second;
170  }
171  tables[this->cur_quad].clear();
172 
173  this->sub_tables = &tables[this->cur_quad];
174  this->update_nodes_ptr();
175 
176  this->order = 20; // fixme
177  }
178 
179  template<typename Scalar>
180  void Filter<Scalar>::free()
181  {
182  for (int i = 0; i < 4; i++)
183  {
184  for(typename std::map<uint64_t, LightArray<struct Filter<Scalar>::Node*>*>::iterator it = tables[i].begin(); it != tables[i].end(); it++)
185  {
186  for(unsigned int l = 0; l < it->second->get_size(); l++)
187  if(it->second->present(l))
188  ::free(it->second->get(l));
189  delete it->second;
190  }
191  tables[i].clear();
192  }
193 
194  if(unimesh)
195  {
196  delete this->mesh;
197  for (int i = 0; i < num; i++)
198  ::free(unidata[i]);
199  delete [] unidata;
200  }
201  }
202 
203  template<typename Scalar>
204  void Filter<Scalar>::reinit()
205  {
206  free();
207  init();
208  }
209 
210  template<typename Scalar>
211  void Filter<Scalar>::push_transform(int son)
212  {
213  MeshFunction<Scalar>::push_transform(son);
214  for (int i = 0; i < num; i++)
215  {
216  // sln_sub[i] contains the value sln[i]->sub_idx, which the Filter thinks
217  // the solution has, or at least had last time. If the filter graph is
218  // cyclic, it could happen that some solutions would get pushed the transform
219  // more than once. This mechanism prevents it. If the filter sees that the
220  // solution already has a different sub_idx than it thinks it should have,
221  // it assumes someone else has already pushed the correct transform. This
222  // is actually the case for cyclic filter graphs and filters used in multi-mesh
223  // computation.
224 
225  if(this->sln[i]->get_transform() == sln_sub[i])
226  this->sln[i]->push_transform(son);
227  sln_sub[i] = this->sln[i]->get_transform();
228  }
229  }
230 
231  template<typename Scalar>
232  void Filter<Scalar>::pop_transform()
233  {
234  MeshFunction<Scalar>::pop_transform();
235  for (int i = 0; i < num; i++)
236  {
237  if(this->sln[i]->get_transform() == sln_sub[i] && sln_sub[i] != 0)
238  this->sln[i]->pop_transform();
239  sln_sub[i] = this->sln[i]->get_transform();
240  }
241  }
242 
243  template<typename Scalar>
244  SimpleFilter<Scalar>::SimpleFilter() : Filter<Scalar>()
245  {
246  }
247 
248  template<typename Scalar>
249  SimpleFilter<Scalar>::SimpleFilter(const Hermes::vector<MeshFunction<Scalar>*>& solutions, const Hermes::vector<int>& items)
250  {
251  this->num = solutions.size();
252  if(this->num > 10)
253  throw Hermes::Exceptions::Exception("Attempt to create an instance of Filter with more than 10 MeshFunctions.");
254  if(items.size() != (unsigned) this->num)
255  if(items.size() > 0)
256  throw Hermes::Exceptions::Exception("Attempt to create an instance of SimpleFilter with different supplied number of MeshFunctions than the number of types of data used from them.");
257 
258  for(int i = 0; i < this->num; i++)
259 
260  {
261  this->sln[i] = solutions.at(i);
262  if(items.size() > 0)
263  this->item[i] = items.at(i);
264  else
265  this->item[i] = H2D_FN_VAL;
266  }
267  this->init();
268  init_components();
269  this->deleteSolutions = false;
270  }
271 
272  template<typename Scalar>
273  SimpleFilter<Scalar>::SimpleFilter(const Hermes::vector<Solution<Scalar>*>& solutions, const Hermes::vector<int>& items)
274  {
275  this->num = solutions.size();
276  if(this->num > 10)
277  throw Hermes::Exceptions::Exception("Attempt to create an instance of Filter with more than 10 MeshFunctions.");
278  if(items.size() != (unsigned) this->num)
279  if(items.size() > 0)
280  throw Hermes::Exceptions::Exception("Attempt to create an instance of SimpleFilter with different supplied number of MeshFunctions than the number of types of data used from them.");
281 
282  for(int i = 0; i < this->num; i++)
283  {
284  this->sln[i] = solutions.at(i);
285  if(items.size() > 0)
286  this->item[i] = items.at(i);
287  else
288  this->item[i] = H2D_FN_VAL;
289  }
290  this->init();
291  init_components();
292  this->deleteSolutions = false;
293  }
294 
295  template<typename Scalar>
296  SimpleFilter<Scalar>::~SimpleFilter()
297  {
298  }
299 
300  template<typename Scalar>
301  void SimpleFilter<Scalar>::init_components()
302  {
303  bool vec1 = false, vec2 = false;
304  for (int i = 0; i < this->num; i++)
305  {
306  if(this->sln[i]->get_num_components() > 1) vec1 = true;
307  if((item[i] & H2D_FN_COMPONENT_0) && (item[i] & H2D_FN_COMPONENT_1)) vec2 = true;
308  if(this->sln[i]->get_num_components() == 1) item[i] &= H2D_FN_COMPONENT_0;
309  }
310  this->num_components = (vec1 && vec2) ? 2 : 1;
311  }
312 
313  template<typename Scalar>
314  void SimpleFilter<Scalar>::precalculate(int order, int mask)
315  {
316  if(mask & (H2D_FN_DX | H2D_FN_DY | H2D_FN_DXX | H2D_FN_DYY | H2D_FN_DXY))
317  throw Hermes::Exceptions::Exception("Filter not defined for derivatives.");
318 
319  Quad2D* quad = this->quads[this->cur_quad];
320  int np = quad->get_num_points(order, this->element->get_mode());
321  struct Function<Scalar>::Node* node = this->new_node(H2D_FN_VAL, np);
322 
323  // precalculate all solutions
324  for (int i = 0; i < this->num; i++)
325  this->sln[i]->set_quad_order(order, item[i]);
326 
327  for (int j = 0; j < this->num_components; j++)
328  {
329  // obtain corresponding tables
330  Scalar* tab[H2D_MAX_COMPONENTS];
331  for (int i = 0; i < this->num; i++)
332  {
333  int a = 0, b = 0, mask = item[i];
334  if(mask >= 0x40) { a = 1; mask >>= 6; }
335  while (!(mask & 1)) { mask >>= 1; b++; }
336  tab[i] = this->sln[i]->get_values(this->num_components == 1 ? a : j, b);
337  if(tab[i] == NULL) throw Hermes::Exceptions::Exception("Value of 'item%d' is incorrect in filter definition.", i + 1);
338  }
339 
340  Hermes::vector<Scalar*> values;
341  for(int i = 0; i < this->num; i++)
342  values.push_back(tab[i]);
343 
344  // apply the filter
345  filter_fn(np, values, node->values[j][0]);
346  }
347 
348  if(this->nodes->present(order))
349  {
350  assert(this->nodes->get(order) == this->cur_node);
351  ::free(this->nodes->get(order));
352  }
353  this->nodes->add(node, order);
354  this->cur_node = node;
355  }
356 
357  template<typename Scalar>
358  Func<Scalar>* SimpleFilter<Scalar>::get_pt_value(double x, double y)
359  {
360  Scalar val[H2D_MAX_COMPONENTS];
361  for (int i = 0; i < this->num; i++)
362  val[i] = this->sln[i]->get_pt_value(x, y)->val[0];
363 
364  Func<Scalar>* toReturn = new Func<Scalar>(1, 1);
365 
366  Scalar result;
367  Hermes::vector<Scalar*> values;
368  for(int i = 0; i < this->num; i++)
369  values.push_back(&val[i]);
370 
371  // apply the filter
372  filter_fn(1, values, &result);
373 
374  toReturn->val[0] = result;
375  return toReturn;
376  }
377 
378  template<typename Scalar>
379  DXFilter<Scalar>::DXFilter(const Hermes::vector<MeshFunction<Scalar>*>& solutions) : DXDYFilter<Scalar>(solutions)
380  {
381  }
382 
383  template<typename Scalar>
384  DXFilter<Scalar>::~DXFilter()
385  {
386  }
387 
388  template<typename Scalar>
389  void DXFilter<Scalar>::filter_fn(int n, Hermes::vector<Scalar *> values, Hermes::vector<Scalar *> dx, Hermes::vector<Scalar *> dy, Scalar* rslt, Scalar* rslt_dx, Scalar* rslt_dy)
390  {
391  for (int i = 0; i < n; i++)
392  {
393  rslt_dx[i] = 0;
394  rslt_dy[i] = 0;
395  for(unsigned int j = 0; j < values.size(); j++)
396  {
397  rslt_dx[i] += dx.at(j)[i];
398  rslt_dy[i] += dy.at(j)[i];
399  }
400  }
401  }
402 
403  template<typename Scalar>
404  MeshFunction<Scalar>* DXFilter<Scalar>::clone() const
405  {
406  Hermes::vector<MeshFunction<Scalar>*> slns;
407  Hermes::vector<int> items;
408  for(int i = 0; i < this->num; i++)
409  {
410  slns.push_back(this->sln[i]->clone());
411  }
412  DXFilter<Scalar>* filter = new DXFilter<Scalar>(slns);
413  filter->setDeleteSolutions();
414  return filter;
415  }
416 
417  ComplexFilter::ComplexFilter(MeshFunction<std::complex<double> >* solution, int item) : Filter<double>()
418  {
419  this->num = 0;
420  this->unimesh = false;
421  this->sln_complex = solution;
422  this->num_components = solution->get_num_components();
423  this->mesh = solution->get_mesh();
424  set_quad_2d(&g_quad_2d_std);
425 
426  // Set NOT to delete solution, as it is probably taken care of by the user.
427  this->deleteSolutions = false;
428  }
429 
430  ComplexFilter::~ComplexFilter()
431  {
432  this->free();
433  }
434 
435  void ComplexFilter::free()
436  {
437  for (int i = 0; i < 4; i++)
438  {
439 #ifdef _MSC_VER // For Visual Studio compiler the latter does not compile.
440  for(std::map<uint64_t, LightArray<Node*>*>::iterator it = tables[i].begin(); it != tables[i].end(); it++)
441 #else
442  for(typename std::map<uint64_t, LightArray<struct Function<double>::Node*>*>::iterator it = tables[i].begin(); it != tables[i].end(); it++)
443 #endif
444  {
445  for(unsigned int l = 0; l < it->second->get_size(); l++)
446  if(it->second->present(l))
447  ::free(it->second->get(l));
448  delete it->second;
449  }
450  tables[i].clear();
451  }
452 
453  if(this->deleteSolutions)
454  delete this->sln_complex;
455  }
456 
457  void ComplexFilter::set_quad_2d(Quad2D* quad_2d)
458  {
459  MeshFunction<double>::set_quad_2d(quad_2d);
460  this->sln_complex->set_quad_2d(quad_2d);
461  }
462 
463  void ComplexFilter::set_active_element(Element* e)
464  {
465  MeshFunction<double>::set_active_element(e);
466 
467  this->sln_complex->set_active_element(e);
468 
469  memset(sln_sub, 0, sizeof(sln_sub));
470 
471  for(std::map<uint64_t, LightArray<struct Function<double>::Node*>*>::iterator it = tables[this->cur_quad].begin(); it != tables[this->cur_quad].end(); it++)
472  {
473  for(unsigned int l = 0; l < it->second->get_size(); l++)
474  if(it->second->present(l))
475  ::free(it->second->get(l));
476  delete it->second;
477  }
478  tables[this->cur_quad].clear();
479 
480  this->sub_tables = &tables[this->cur_quad];
481 
482  this->update_nodes_ptr();
483 
484  this->order = 20; // fixme
485  }
486 
487  void ComplexFilter::push_transform(int son)
488  {
489  MeshFunction<double>::push_transform(son);
490  this->sln_complex->push_transform(son);
491  }
492 
493  void ComplexFilter::pop_transform()
494  {
495  MeshFunction<double>::pop_transform();
496  this->sln_complex->pop_transform();
497  }
498 
499  void ComplexFilter::precalculate(int order, int mask)
500  {
501  if(mask & (H2D_FN_DX | H2D_FN_DY | H2D_FN_DXX | H2D_FN_DYY | H2D_FN_DXY))
502  throw Hermes::Exceptions::Exception("Filter not defined for derivatives.");
503 
504  Quad2D* quad = this->quads[this->cur_quad];
505  int np = quad->get_num_points(order, this->element->get_mode());
506  struct Function<double>::Node* node = this->new_node(H2D_FN_VAL, np);
507 
508  this->sln_complex->set_quad_order(order, H2D_FN_VAL);
509 
510  // obtain corresponding tables
511  filter_fn(np, this->sln_complex->get_values(0, 0), node->values[0][0]);
512  if(num_components > 1)
513  filter_fn(np, this->sln_complex->get_values(1, 0), node->values[1][0]);
514 
515  if(this->nodes->present(order))
516  {
517  assert(this->nodes->get(order) == this->cur_node);
518  ::free(this->nodes->get(order));
519  }
520 
521  this->cur_node = node;
522  }
523 
524  Func<double>* ComplexFilter::get_pt_value(double x, double y)
525  {
526  Func<std::complex<double> >* val = this->sln_complex->get_pt_value(x, y);
527 
528  Func<double>* toReturn = new Func<double>(1, this->num_components);
529 
530  double result;
531 
532  // apply the filter
533  filter_fn(1, val->val, &result);
534 
535  if(this->num_components == 1)
536  {
537  toReturn->val[0] = result;
538  toReturn->dx[0] = 0.0;
539  toReturn->dy[0] = 0.0;
540  }
541  else
542  {
543  this->warn("ComplexFilter only implemented for scalar functions.");
544  }
545  return toReturn;
546  }
547 
548  template<typename Scalar>
549  DXDYFilter<Scalar>::DXDYFilter() : Filter<Scalar>()
550  {
551  }
552 
553  template<typename Scalar>
554  DXDYFilter<Scalar>::DXDYFilter(const Hermes::vector<MeshFunction<Scalar>*>& solutions) : Filter<Scalar>(solutions)
555  {
556  init_components();
557  }
558 
559  template<typename Scalar>
560  DXDYFilter<Scalar>::DXDYFilter(const Hermes::vector<Solution<Scalar>*>& solutions) : Filter<Scalar>(solutions)
561  {
562  init_components();
563  }
564 
565  template<typename Scalar>
566  DXDYFilter<Scalar>::~DXDYFilter()
567  {
568  }
569 
570  template<typename Scalar>
571  void DXDYFilter<Scalar>::init_components()
572  {
573  this->num_components = this->sln[0]->get_num_components();
574  for (int i = 1; i < this->num; i++)
575  if(this->sln[i]->get_num_components() != this->num_components)
576  throw Hermes::Exceptions::Exception("Filter: Solutions do not have the same number of components!");
577  }
578 
579  template<typename Scalar>
580  void DXDYFilter<Scalar>::init(const Hermes::vector<MeshFunction<Scalar>*>& solutions)
581  {
582  Filter<Scalar>::init(solutions);
583  init_components();
584  }
585 
586  template<typename Scalar>
587  void DXDYFilter<Scalar>::precalculate(int order, int mask)
588  {
589  Quad2D* quad = this->quads[this->cur_quad];
590  int np = quad->get_num_points(order, this->element->get_mode());
591  struct Function<Scalar>::Node* node = this->new_node(H2D_FN_DEFAULT, np);
592 
593  // precalculate all solutions
594  for (int i = 0; i < this->num; i++)
595  this->sln[i]->set_quad_order(order, H2D_FN_DEFAULT);
596 
597  for (int j = 0; j < this->num_components; j++)
598  {
599  // obtain solution tables
600  Scalar *val[H2D_MAX_COMPONENTS], *dx[H2D_MAX_COMPONENTS], *dy[H2D_MAX_COMPONENTS];
601  for (int i = 0; i < this->num; i++)
602  {
603  val[i] = this->sln[i]->get_fn_values(j);
604  dx[i] = this->sln[i]->get_dx_values(j);
605  dy[i] = this->sln[i]->get_dy_values(j);
606  }
607 
608  Hermes::vector<Scalar *> values_vector;
609  Hermes::vector<Scalar *> dx_vector;
610  Hermes::vector<Scalar *> dy_vector;
611 
612  for(int i = 0; i < this->num; i++)
613  {
614  values_vector.push_back(val[i]);
615  dx_vector.push_back(dx[i]);
616  dy_vector.push_back(dy[i]);
617  }
618 
619  // apply the filter
620  filter_fn(np, values_vector, dx_vector, dy_vector, node->values[j][0], node->values[j][1], node->values[j][2]);
621  }
622 
623  if(this->nodes->present(order))
624  {
625  assert(this->nodes->get(order) == this->cur_node);
626  ::free(this->nodes->get(order));
627  }
628  this->nodes->add(node, order);
629  this->cur_node = node;
630  }
631 
632  template<typename Scalar>
633  Func<Scalar>* DXDYFilter<Scalar>::get_pt_value(double x, double y)
634  {
635  return 0;
636  }
637 
638  template<typename Scalar>
639  void MagFilter<Scalar>::filter_fn(int n, Hermes::vector<Scalar*> values, Scalar* result)
640  {
641  for (int i = 0; i < n; i++)
642 
643  {
644  result[i] = 0;
645  for(unsigned int j = 0; j < values.size(); j++)
646  result[i] += sqr(values.at(j)[i]);
647  result[i] = sqrt(result[i]);
648  }
649  };
650 
651  template<typename Scalar>
652  MagFilter<Scalar>::MagFilter(Hermes::vector<MeshFunction<Scalar>*> solutions, Hermes::vector<int> items) : SimpleFilter<Scalar>(solutions, items)
653  {
654  };
655 
656  template<typename Scalar>
657  MagFilter<Scalar>::MagFilter(MeshFunction<Scalar>* sln1, int item1)
658  : SimpleFilter<Scalar>(Hermes::vector<MeshFunction<Scalar>*>(sln1, sln1),
659  Hermes::vector<int>(item1 & H2D_FN_COMPONENT_0, item1 & H2D_FN_COMPONENT_1))
660  {
661  if(sln1->get_num_components() < 2)
662  throw Hermes::Exceptions::Exception("The single-argument constructor is intended for vector-valued solutions.");
663  };
664 
665  template<typename Scalar>
666  MagFilter<Scalar>::~MagFilter()
667  {
668  };
669 
670  template<typename Scalar>
671  MeshFunction<Scalar>* MagFilter<Scalar>::clone() const
672  {
673  Hermes::vector<MeshFunction<Scalar>*> slns;
674  Hermes::vector<int> items;
675  for(int i = 0; i < this->num; i++)
676  {
677  slns.push_back(this->sln[i]->clone());
678  items.push_back(this->item[i]);
679  }
680  MagFilter<Scalar>* filter = new MagFilter<Scalar>(slns, items);
681  filter->setDeleteSolutions();
682  return filter;
683  }
684 
685  void TopValFilter::filter_fn(int n, Hermes::vector<double*> values, double* result)
686  {
687  for (int i = 0; i < n; i++)
688  {
689  result[i] = 0;
690  for(unsigned int j = 0; j < values.size(); j++)
691  if(values.at(j)[i] > limits[j])
692  result[i] = limits[j];
693  else
694  result[i] = values.at(j)[i];
695  }
696  };
697 
698  TopValFilter::TopValFilter(Hermes::vector<MeshFunction<double>*> solutions, Hermes::vector<double> limits, Hermes::vector<int> items) : SimpleFilter<double>(solutions, items), limits(limits)
699  {
700  };
701 
702  TopValFilter::TopValFilter(MeshFunction<double>* sln, double limit, int item)
703  : SimpleFilter<double>()
704  {
705  this->limits.push_back(limit);
706  this->sln[0] = sln;
707  this->item[0] = item;
708  this->num = 1;
709  Filter<double>::init();
710  };
711 
712  TopValFilter::~TopValFilter()
713  {
714  }
715 
716  MeshFunction<double>* TopValFilter::clone() const
717  {
718  Hermes::vector<MeshFunction<double>*> slns;
719  Hermes::vector<int> items;
720  for(int i = 0; i < this->num; i++)
721  {
722  slns.push_back(this->sln[i]->clone());
723  items.push_back(this->item[i]);
724  }
725  TopValFilter* filter = new TopValFilter(slns, limits, items);
726  filter->setDeleteSolutions();
727  return filter;
728  }
729 
730  void BottomValFilter::filter_fn(int n, Hermes::vector<double*> values, double* result)
731  {
732  for (int i = 0; i < n; i++)
733  {
734  result[i] = 0;
735  for(unsigned int j = 0; j < values.size(); j++)
736  if(values.at(j)[i] < limits[j])
737  result[i] = limits[j];
738  else
739  result[i] = values.at(j)[i];
740  }
741  };
742 
743  BottomValFilter::BottomValFilter(Hermes::vector<MeshFunction<double>*> solutions, Hermes::vector<double> limits, Hermes::vector<int> items) : SimpleFilter<double>(solutions, items), limits(limits)
744  {
745  };
746 
747  BottomValFilter::BottomValFilter(MeshFunction<double>* sln, double limit, int item)
748  : SimpleFilter<double>()
749  {
750  this->limits.push_back(limit);
751  this->sln[0] = sln;
752  this->item[0] = item;
753  this->num = 1;
754  Filter<double>::init();
755  };
756 
757  BottomValFilter::~BottomValFilter()
758  {
759  }
760 
761  MeshFunction<double>* BottomValFilter::clone() const
762  {
763  Hermes::vector<MeshFunction<double>*> slns;
764  Hermes::vector<int> items;
765  for(int i = 0; i < this->num; i++)
766  {
767  slns.push_back(this->sln[i]->clone());
768  items.push_back(this->item[i]);
769  }
770  BottomValFilter* filter = new BottomValFilter(slns, limits, items);
771  filter->setDeleteSolutions();
772  return filter;
773  }
774 
775  void ValFilter::filter_fn(int n, Hermes::vector<double*> values, double* result)
776  {
777  for (int i = 0; i < n; i++)
778  {
779  result[i] = 0;
780  for(unsigned int j = 0; j < values.size(); j++)
781  if(values.at(j)[i] < low_limits[j])
782  result[i] = low_limits[j];
783  else
784  if(values.at(j)[i] > high_limits[j])
785  result[i] = high_limits[j];
786  else
787  result[i] = values.at(j)[i];
788  }
789  };
790 
791  ValFilter::ValFilter(Hermes::vector<MeshFunction<double>*> solutions, Hermes::vector<double> low_limits, Hermes::vector<double> high_limits, Hermes::vector<int> items) : SimpleFilter<double>(solutions, items), low_limits(low_limits), high_limits(high_limits)
792  {
793  };
794 
795  ValFilter::ValFilter(MeshFunction<double>* sln, double low_limit, double high_limit, int item)
796  : SimpleFilter<double>()
797  {
798  this->low_limits.push_back(low_limit);
799  this->high_limits.push_back(high_limit);
800  this->sln[0] = sln;
801  this->item[0] = item;
802  this->num = 1;
803  Filter<double>::init();
804  };
805 
806  ValFilter::~ValFilter()
807  {
808  }
809 
810  MeshFunction<double>* ValFilter::clone() const
811  {
812  Hermes::vector<MeshFunction<double>*> slns;
813  Hermes::vector<int> items;
814  for(int i = 0; i < this->num; i++)
815  {
816  slns.push_back(this->sln[i]->clone());
817  items.push_back(this->item[i]);
818  }
819  ValFilter* filter = new ValFilter(slns, low_limits, high_limits, items);
820  filter->setDeleteSolutions();
821  return filter;
822  }
823 
824  template<typename Scalar>
825  void DiffFilter<Scalar>::filter_fn(int n, Hermes::vector<Scalar*> values, Scalar* result)
826  {
827  for (int i = 0; i < n; i++) result[i] = values.at(0)[i] - values.at(1)[i];
828  };
829 
830  template<typename Scalar>
831  DiffFilter<Scalar>::DiffFilter(Hermes::vector<MeshFunction<Scalar>*> solutions, Hermes::vector<int> items) : SimpleFilter<Scalar>(solutions, items) {}
832 
833  template<typename Scalar>
834  DiffFilter<Scalar>::~DiffFilter()
835  {
836  }
837 
838  template<typename Scalar>
839  MeshFunction<Scalar>* DiffFilter<Scalar>::clone() const
840  {
841  Hermes::vector<MeshFunction<Scalar>*> slns;
842  Hermes::vector<int> items;
843  for(int i = 0; i < this->num; i++)
844  {
845  slns.push_back(this->sln[i]->clone());
846  items.push_back(this->item[i]);
847  }
848  DiffFilter* filter = new DiffFilter<Scalar>(slns, items);
849  filter->setDeleteSolutions();
850  return filter;
851  }
852 
853  template<typename Scalar>
854  void SumFilter<Scalar>::filter_fn(int n, Hermes::vector<Scalar*> values, Scalar* result)
855  {
856  for (int i = 0; i < n; i++)
857  {
858  result[i] = 0;
859  for (unsigned int j = 0; j < values.size(); j++)
860  result[i] += values.at(j)[i];
861  }
862  };
863 
864  template<typename Scalar>
865  SumFilter<Scalar>::SumFilter(Hermes::vector<MeshFunction<Scalar>*> solutions, Hermes::vector<int> items) : SimpleFilter<Scalar>(solutions, items) {}
866 
867  template<typename Scalar>
868  SumFilter<Scalar>::~SumFilter()
869  {
870  }
871 
872  template<typename Scalar>
873  MeshFunction<Scalar>* SumFilter<Scalar>::clone() const
874  {
875  Hermes::vector<MeshFunction<Scalar>*> slns;
876  Hermes::vector<int> items;
877  for(int i = 0; i < this->num; i++)
878  {
879  slns.push_back(this->sln[i]->clone());
880  items.push_back(this->item[i]);
881  }
882  SumFilter<Scalar>* filter = new SumFilter<Scalar>(slns, items);
883  filter->setDeleteSolutions();
884  return filter;
885  }
886 
887  template<>
888  void SquareFilter<double>::filter_fn(int n, Hermes::vector<double *> v1, double* result)
889  {
890  for (int i = 0; i < n; i++)
891  result[i] = sqr(v1.at(0)[i]);
892  };
893 
894  template<>
895  void SquareFilter<std::complex<double> >::filter_fn(int n, Hermes::vector<std::complex<double> *> v1, std::complex<double> * result)
896  {
897  for (int i = 0; i < n; i++)
898  result[i] = std::norm(v1.at(0)[i]);
899  };
900 
901  template<typename Scalar>
902  SquareFilter<Scalar>::SquareFilter(Hermes::vector<MeshFunction<Scalar>*> solutions, Hermes::vector<int> items)
903  : SimpleFilter<Scalar>(solutions, items)
904  {
905  if(solutions.size() > 1)
906  throw Hermes::Exceptions::Exception("SquareFilter only supports one MeshFunction.");
907  };
908 
909  template<typename Scalar>
910  SquareFilter<Scalar>::~SquareFilter()
911  {
912  }
913 
914  template<typename Scalar>
915  MeshFunction<Scalar>* SquareFilter<Scalar>::clone() const
916  {
917  Hermes::vector<MeshFunction<Scalar>*> slns;
918  Hermes::vector<int> items;
919  for(int i = 0; i < this->num; i++)
920  {
921  slns.push_back(this->sln[i]->clone());
922  items.push_back(this->item[i]);
923  }
924  SquareFilter<Scalar>* filter = new SquareFilter<Scalar>(slns, items);
925  filter->setDeleteSolutions();
926  return filter;
927  }
928 
929  void AbsFilter::filter_fn(int n, Hermes::vector<double*> v1, double * result)
930  {
931  for (int i = 0; i < n; i++)
932  result[i] = std::abs(v1.at(0)[i]);
933  };
934 
935  AbsFilter::AbsFilter(Hermes::vector<MeshFunction<double>*> solutions, Hermes::vector<int> items)
936  : SimpleFilter<double>(solutions, items)
937  {
938  if(solutions.size() > 1)
939  throw Hermes::Exceptions::Exception("AbsFilter only supports one MeshFunction.");
940  };
941 
942  AbsFilter::AbsFilter(MeshFunction<double>* solution)
943  : SimpleFilter<double>()
944  {
945  this->num = 1;
946  this->sln[0] = solution;
947 
948  this->item[0] = H2D_FN_VAL;
949 
950  this->init();
951  init_components();
952  };
953 
954  AbsFilter::~AbsFilter()
955  {
956  }
957 
958  MeshFunction<double>* AbsFilter::clone() const
959  {
960  Hermes::vector<MeshFunction<double>*> slns;
961  Hermes::vector<int> items;
962  for(int i = 0; i < this->num; i++)
963  {
964  slns.push_back(this->sln[i]->clone());
965  items.push_back(this->item[i]);
966  }
967  AbsFilter* filter = new AbsFilter(slns, items);
968  filter->setDeleteSolutions();
969  return filter;
970  }
971 
972  void RealFilter::filter_fn(int n, std::complex<double>* values, double* result)
973  {
974  for (int i = 0; i < n; i++)
975  result[i] = values[i].real();
976  };
977 
978  MeshFunction<double>* RealFilter::clone() const
979  {
980  RealFilter* filter = new RealFilter(this->sln_complex->clone(), this->item);
981  filter->setDeleteSolutions();
982  return filter;
983  }
984 
985  RealFilter::RealFilter(MeshFunction<std::complex<double> >* solution, int item)
986  : ComplexFilter(solution, item)
987  {
988  };
989 
990  RealFilter::~RealFilter()
991  {
992  }
993 
994  void ImagFilter::filter_fn(int n, std::complex<double>* values, double* result)
995  {
996  for (int i = 0; i < n; i++)
997  result[i] = values[i].imag();
998  };
999 
1000  ImagFilter::ImagFilter(MeshFunction<std::complex<double> >* solution, int item)
1001  : ComplexFilter(solution, item)
1002  {
1003  };
1004 
1005  ImagFilter::~ImagFilter()
1006  {
1007  }
1008 
1009  MeshFunction<double>* ImagFilter::clone() const
1010  {
1011  ImagFilter* filter = new ImagFilter(this->sln_complex->clone(), this->item);
1012  filter->setDeleteSolutions();
1013  return filter;
1014  }
1015 
1016  void ComplexAbsFilter::filter_fn(int n, std::complex<double>* values, double* result)
1017  {
1018  for (int i = 0; i < n; i++)
1019  result[i] = sqrt(sqr(values[i].real()) + sqr(values[i].imag()));
1020  };
1021 
1022  MeshFunction<double>* ComplexAbsFilter::clone() const
1023  {
1024  ComplexAbsFilter* filter = new ComplexAbsFilter(this->sln_complex->clone(), this->item);
1025  filter->setDeleteSolutions();
1026  return filter;
1027  }
1028 
1029  ComplexAbsFilter::ComplexAbsFilter(MeshFunction<std::complex<double> >* solution, int item)
1030  : ComplexFilter(solution, item)
1031  {
1032  };
1033 
1034  ComplexAbsFilter::~ComplexAbsFilter()
1035  {
1036  }
1037 
1038  void AngleFilter::filter_fn(int n, Hermes::vector<std::complex<double>*> v1, double* result)
1039  {
1040  for (int i = 0; i < n; i++)
1041  result[i] = atan2( v1.at(0)[i].imag(), v1.at(0)[i].real() );
1042  };
1043 
1044  AngleFilter::AngleFilter(Hermes::vector<MeshFunction<std::complex<double> >*> solutions, Hermes::vector<int> items)
1045  : SimpleFilter<std::complex<double> >(solutions, items)
1046  {
1047  if(solutions.size() > 1)
1048  throw Hermes::Exceptions::Exception("RealFilter only supports one MeshFunction.");
1049  };
1050 
1051  AngleFilter::~AngleFilter()
1052  {
1053  }
1054 
1055  void VonMisesFilter::precalculate(int order, int mask)
1056  {
1057  if(mask & (H2D_FN_DX | H2D_FN_DY | H2D_FN_DXX | H2D_FN_DYY | H2D_FN_DXY))
1058  throw Hermes::Exceptions::Exception("VonMisesFilter not defined for derivatives.");
1059 
1060  Quad2D* quad = this->quads[this->cur_quad];
1061  int np = quad->get_num_points(order, this->element->get_mode());
1062  Filter<double>::Node* node = new_node(H2D_FN_VAL_0, np);
1063 
1064  this->sln[0]->set_quad_order(order, H2D_FN_VAL | H2D_FN_DX | H2D_FN_DY);
1065  this->sln[1]->set_quad_order(order, H2D_FN_DX | H2D_FN_DY);
1066 
1067  double *dudx, *dudy, *dvdx, *dvdy;
1068  this->sln[0]->get_dx_dy_values(dudx, dudy);
1069  this->sln[1]->get_dx_dy_values(dvdx, dvdy);
1070  double *uval = this->sln[0]->get_fn_values();
1071  update_refmap();
1072  double *x = refmap->get_phys_x(order);
1073 
1074  for (int i = 0; i < np; i++)
1075  {
1076  // stress tensor
1077  double tz = lambda*(dudx[i] + dvdy[i]);
1078  double tx = tz + 2*mu*dudx[i];
1079  double ty = tz + 2*mu*dvdy[i];
1080  if(cyl) tz += 2*mu*uval[i] / x[i];
1081  double txy = mu*(dudy[i] + dvdx[i]);
1082 
1083  // Von Mises stress
1084  node->values[0][0][i] = 1.0/sqrt(2.0) * sqrt(sqr(tx - ty) + sqr(ty - tz) + sqr(tz - tx) + 6*sqr(txy));
1085  }
1086 
1087  if(this->nodes->present(order))
1088  {
1089  assert(this->nodes->get(order) == cur_node);
1090  ::free(this->nodes->get(order));
1091  }
1092  this->nodes->add(node, order);
1093  cur_node = node;
1094  }
1095 
1096  Func<double>* VonMisesFilter::get_pt_value(double x, double y)
1097  {
1098  return 0;
1099  }
1100 
1101  VonMisesFilter::VonMisesFilter(Hermes::vector<MeshFunction<double>*> solutions, double lambda, double mu,
1102  int cyl, int item1, int item2)
1103  : Filter<double>(solutions)
1104  {
1105  this->mu = mu;
1106  this->lambda = lambda;
1107  this->cyl = cyl;
1108  this->item1 = item1;
1109  this->item2 = item2;
1110  }
1111 
1112  VonMisesFilter::VonMisesFilter(MeshFunction<double>** solutions, int num, double lambda, double mu,
1113  int cyl, int item1, int item2): Filter<double>(solutions, num)
1114  {
1115  this->mu = mu;
1116  this->lambda = lambda;
1117  this->cyl = cyl;
1118  this->item1 = item1;
1119  this->item2 = item2;
1120  }
1121 
1122  VonMisesFilter::~VonMisesFilter()
1123  {
1124  }
1125 
1126  MeshFunction<double>* VonMisesFilter::clone() const
1127  {
1128  MeshFunction<double>** slns = new MeshFunction<double>*[num];
1129  for(int i = 0; i < num; i++)
1130  slns[i] = sln[i]->clone();
1131  VonMisesFilter* filter = new VonMisesFilter(slns, num, lambda, mu, cyl, item1, item2);
1132  filter->setDeleteSolutions();
1133  return filter;
1134  }
1135 
1136  template<typename Scalar>
1137  void LinearFilter<Scalar>::precalculate(int order, int mask)
1138  {
1139  Quad2D* quad = this->quads[this->cur_quad];
1140  int np = quad->get_num_points(order, this->element->get_mode());
1141  struct Filter<Scalar>::Node* node = this->new_node(H2D_FN_DEFAULT, np);
1142 
1143  // precalculate all solutions
1144  for (int i = 0; i < this->num; i++)
1145  this->sln[i]->set_quad_order(order);
1146 
1147  for (int j = 0; j < this->num_components; j++)
1148  {
1149  // obtain solution tables
1150  Scalar *val[H2D_MAX_COMPONENTS], *dx[H2D_MAX_COMPONENTS], *dy[H2D_MAX_COMPONENTS];
1151  for (int i = 0; i < this->num; i++)
1152  {
1153  val[i] = this->sln[i]->get_fn_values(j);
1154  dx[i] = this->sln[i]->get_dx_values(j);
1155  dy[i] = this->sln[i]->get_dy_values(j);
1156  }
1157  if(this->num == 2)
1158  for (int i = 0; i < np; i++)
1159  {
1160  node->values[j][0][i] = tau_frac * (val[1][i] - val[0][i]) + val[1][i];
1161  node->values[j][1][i] = tau_frac * (dx[1][i] - dx[0][i]) + dx[1][i];
1162  node->values[j][2][i] = tau_frac * (dy[1][i] - dy[0][i]) + dy[1][i];
1163  }
1164  else
1165  for (int i = 0; i < np; i++)
1166  {
1167  node->values[j][0][i] = val[0][i];
1168  node->values[j][1][i] = dx[0][i];
1169  node->values[j][2][i] = dy[0][i];
1170  }
1171  }
1172 
1173  if(this->nodes->present(order))
1174  {
1175  assert(this->nodes->get(order) == this->cur_node);
1176  ::free(this->nodes->get(order));
1177  }
1178  this->nodes->add(node, order);
1179  this->cur_node = node;
1180  }
1181 
1182  template<typename Scalar>
1183  Func<Scalar>* LinearFilter<Scalar>::get_pt_value(double x, double y)
1184  {
1185  return 0;
1186  }
1187 
1188  template<typename Scalar>
1189  LinearFilter<Scalar>::LinearFilter(MeshFunction<Scalar>* old) // : Filter(old)
1190  {
1191  init_components();
1192  }
1193 
1194  template<typename Scalar>
1195  LinearFilter<Scalar>::LinearFilter(MeshFunction<Scalar>* older, MeshFunction<Scalar>* old, double tau_frac)
1196  : Filter<Scalar>(Hermes::vector<MeshFunction<Scalar>*>(older, old))
1197  {
1198  this->tau_frac = tau_frac;
1199  init_components();
1200  }
1201 
1202  template<typename Scalar>
1203  LinearFilter<Scalar>::~LinearFilter()
1204  {
1205  }
1206 
1207  template<typename Scalar>
1208  void LinearFilter<Scalar>::init_components()
1209  {
1210  this->num_components = this->sln[0]->get_num_components();
1211  for (int i = 1; i < this->num; i++)
1212  if(this->sln[i]->get_num_components() != this->num_components)
1213  throw Hermes::Exceptions::Exception("Filter: Solutions do not have the same number of components!");
1214  }
1215 
1216  template<typename Scalar>
1217  void LinearFilter<Scalar>::set_active_element(Element* e)
1218  {
1219  Filter<Scalar>::set_active_element(e);
1220 
1221  this->order = 0;
1222  for (int i = 0; i < this->num; i++)
1223  {
1224  int o = this->sln[i]->get_fn_order();
1225  if(o > this->order) this->order = o;
1226  }
1227  }
1228 
1229  template class HERMES_API Filter<double>;
1230  template class HERMES_API Filter<std::complex<double> >;
1231  template class HERMES_API SimpleFilter<double>;
1232  template class HERMES_API SimpleFilter<std::complex<double> >;
1233  template class HERMES_API DXDYFilter<double>;
1234  template class HERMES_API DXDYFilter<std::complex<double> >;
1235  template class HERMES_API DXFilter<double>;
1236  template class HERMES_API DXFilter<std::complex<double> >;
1237  template class HERMES_API MagFilter<double>;
1238  template class HERMES_API MagFilter<std::complex<double> >;
1239  template class HERMES_API DiffFilter<double>;
1240  template class HERMES_API DiffFilter<std::complex<double> >;
1241  template class HERMES_API SumFilter<double>;
1242  template class HERMES_API SumFilter<std::complex<double> >;
1243  template class HERMES_API SquareFilter<double>;
1244  template class HERMES_API SquareFilter<std::complex<double> >;
1245  }
1246 }