← Back to branch summary

~gabriel1984sibiu/agros2d/agros2d

« back to all changes in this revision

Viewing changes to hermes/hermes2d/src/function/filter.cpp

  • Committer: Grevutiu Gabriel
  • Date: 2014-11-15 19:05:36 UTC
  • Revision ID: gabriel1984sibiu@gmail.com-20141115190536-1d4q8ez0f8b89ktj
original upstream code

Show diffs side-by-side

added added

removed removed

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