vorticity_smoother.h
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Version 1.0; svn revision $LastChangedRevision: 1294 $
7 //LIC//
8 //LIC// $LastChangedDate: 2017-04-17 14:08:43 +0100 (Mon, 17 Apr 2017) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 //Header file for Navier Stokes elements
31 
32 #ifndef OOMPH_VORTICITY_SMOOTHER_HEADER
33 #define OOMPH_VORTICITY_SMOOTHER_HEADER
34 
35 // Config header generated by autoconfig
36 #ifdef HAVE_CONFIG_H
37 #include <oomph-lib-config.h>
38 #endif
39 
40 
41 //===============================================================
42 // WARNING: THIS IS WORK IN PROGRESS -- ONLY USED IN 2d SO FAR
43 //===============================================================
44 
45 
46 // namespace extension
47 namespace oomph
48 {
49 
50 //===============================================
51 /// Overloaded element that allows projection of
52 /// vorticity.
53 //===============================================
54 template<class ELEMENT>
55 class VorticitySmootherElement : public virtual ELEMENT
56 {
57 
58 public:
59 
60 
61  /// Constructor
63  {
64  // hierher only coded up for 2D at the moment. Check.
65 
66  // Index of smoothed vorticity
68 
69  // Number of values per field: u,v,p,omega,d/dx,d/dy,
70  // d^2/dx^2,d^2/dxdy, d^2/dy^2,
71  // d^3/dx^3, d^3/dx^2dy, d^3/dxdy^2,d^3,dy^3
72  // du/dx, du/dy, dv/dx, dv/dy
73  Number_of_values_per_field=17; // hierher 13;
74 
75  // Pointer to fct that specifies exact vorticity and
76  // derivs (for validation).
78  }
79 
80  /// Typedef for pointer to function that specifies the exact
81  /// vorticity and derivs (for validation)
82  typedef void (*ExactVorticityFctPt)(const Vector<double>& x,
83  double& vort,
84  Vector<double>& dvort_dx,
85  Vector<double>& dvort_dxdy,
86  Vector<double>& dvort_dxdxdy,
87  Vector<double>& dveloc_dx);
88 
89  /// Access function: Pointer to fct that specifies exact vorticity and
90  /// derivs (for validation).
92  {return Exact_vorticity_fct_pt;}
93 
94  /// Access function: Pointer to fct that specifies exact vorticity and
95  /// derivs (for validation). const version
97  {return Exact_vorticity_fct_pt;}
98 
99 
100  /// Index of smoothed vorticity -- followed by derivatives
101  unsigned smoothed_vorticity_index() const
102  {
104  }
105 
106 
107  /// \short Number of values required at local node n. In order to simplify
108  /// matters, we allocate storage for pressure variables at all the nodes
109  /// and then pin those that are not used.
110  unsigned required_nvalue(const unsigned &n) const
111  {return Number_of_values_per_field;}
112 
113 
114  /// \short Number of continuously interpolated values:
115  unsigned ncont_interpolated_values() const
116  {return Number_of_values_per_field;}
117 
118 /// \short Get the function value u in Vector.
119 /// Note: Given the generality of the interface (this function
120 /// is usually called from black-box documentation or interpolation routines),
121 /// the values Vector sets its own size in here.
123  {
124  unsigned t=0;
125  get_interpolated_values(t,s,values);
126  }
127 
128 /// \short Get the function value u in Vector.
129 /// Note: Given the generality of the interface (this function
130 /// is usually called from black-box documentation or interpolation routines),
131 /// the values Vector sets its own size in here.
132  void get_interpolated_values(const unsigned& t, const Vector<double>&s,
133  Vector<double>& values)
134  {
135  unsigned DIM=2;
136 
137  // Set size of Vector
138  values.resize(Number_of_values_per_field);
139 
140  // Initialise
141  for(unsigned i=0;i<DIM+1;i++) {values[i]=0.0;}
142 
143  //Find out how many nodes there are
144  unsigned n_node = this->nnode();
145 
146  // Shape functions
147  Shape psif(n_node);
148  this->shape(s,psif);
149 
150  //Calculate velocities: values[0],...
151  for(unsigned i=0;i<DIM;i++)
152  {
153  //Get the index at which the i-th velocity is stored
154  unsigned u_nodal_index = this->u_index_nst(i);
155  for(unsigned l=0;l<n_node;l++)
156  {
157  values[i] += this->nodal_value(t,l,u_nodal_index)*psif[l];
158  }
159  }
160 
161  //Calculate pressure: values[DIM]
162  //(no history is carried in the pressure)
163  values[DIM] = this->interpolated_p_nst(s);
164 
165 
166  // No need to interpolate these onto the new mesh
167  /* // Vorticity and its derivs */
168  /* for (unsigned i=0;i<6;i++) */
169  /* { */
170  /* for(unsigned l=0;l<n_node;l++) */
171  /* { */
172  /* values[3+i]+=this->nodal_value(t,l,Smoothed_vorticity_index+i)*psif[l]; */
173  /* } */
174  /* } */
175 
176  }
177 
178  /// Pin all smoothed vorticity quantities
180  {
181  unsigned nnod=this->nnode();
182  for (unsigned j=0;j<nnod;j++)
183  {
184  Node* nod_pt=this->node_pt(j);
185  for (unsigned i=Smoothed_vorticity_index;
187  {
188  nod_pt->pin(i);
189  }
190  }
191  }
192 
193 
194  /// \short Output exact veloc, vorticity, derivs and indicator
195  /// based on functions specified by two function pointers
196  void output_analytical_veloc_and_vorticity(std::ostream &outfile,
197  const unsigned &nplot)
198  {
199  //Vector of local coordinates
200  Vector<double> s(2);
201 
202  // Shape functions
203  unsigned n_node = this->nnode();
204  Shape psif(n_node);
205 
206  // Tecplot header info
207  outfile << this->tecplot_zone_string(nplot);
208 
209  // Loop over plot points
210  unsigned num_plot_points=this->nplot_points(nplot);
211  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
212  {
213  // Get local coordinates of plot point
214  this->get_s_plot(iplot,nplot,s);
215 
216  // Coordinates
217  Vector<double> x(2);
218  for(unsigned i=0;i<2;i++)
219  {
220  x[i]=this->interpolated_x(s,i);
221  outfile << x[i] << " ";
222  }
223 
224  // Fake veloc and pressure
225  outfile << "0.0 0.0 0.0 ";
226 
227  // Get vorticity and its derivatives
228  double vort=0.0;
229  Vector<double> dvort_dx(2);
230  Vector<double> dvort_dxdy(3);
231  Vector<double> dvort_dxdxdy(4);
232  Vector<double> dveloc_dx(4);
234  vort,
235  dvort_dx,
236  dvort_dxdy,
237  dvort_dxdxdy,
238  dveloc_dx);
239 
240  // Smoothed vorticity
241  outfile << vort << " ";
242 
243  // Smoothed vorticity derivatives (d/dx, d/dy, d^2/dx^2, d^2/dxdy, d^2/dy^2
244  // d^3/dx^3, d^3/dx^2dy, d^3/dxdy^2, d^3/dy^3,
245  // du/dx, du/dy, dv/dx, dv/dy
246  outfile << dvort_dx[0] << " "
247  << dvort_dx[1] << " "
248  << dvort_dxdy[0] << " "
249  << dvort_dxdy[1] << " "
250  << dvort_dxdy[2] << " "
251  << dvort_dxdxdy[0] << " "
252  << dvort_dxdxdy[1] << " "
253  << dvort_dxdxdy[2] << " "
254  << dvort_dxdxdy[3] << " "
255  << dveloc_dx[0] << " "
256  << dveloc_dx[1] << " "
257  << dveloc_dx[2] << " "
258  << dveloc_dx[3] << " ";
259 
260  outfile << std::endl;
261  }
262 
263  // Write tecplot footer (e.g. FE connectivity lists)
264  this->write_tecplot_zone_footer(outfile,nplot);
265  }
266 
267 
268 
269 
270  /// \short Output veloc, smoothed vorticity and derivatives
271  void output_smoothed_vorticity(std::ostream &outfile,
272  const unsigned &nplot)
273  {
274  //Vector of local coordinates
275  Vector<double> s(2);
276 
277  // Shape functions
278  unsigned n_node = this->nnode();
279  Shape psif(n_node);
280 
281  // Tecplot header info
282  outfile << this->tecplot_zone_string(nplot);
283 
284  // Loop over plot points
285  unsigned num_plot_points=this->nplot_points(nplot);
286  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
287  {
288  // Get local coordinates of plot point
289  this->get_s_plot(iplot,nplot,s);
290 
291  // Get vorticity and its derivatives (reconstructed)
292  double vort=0.0;
293  Vector<double> veloc(2);
294  Vector<double> dvort_dx(2);
295  Vector<double> dvort_dxdy(3);
296  Vector<double> dvort_dxdxdy(4);
297  Vector<double> dveloc_dx(4);
299  veloc,
300  vort,
301  dvort_dx,
302  dvort_dxdy,
303  dvort_dxdxdy,
304  dveloc_dx);
305  // Coordinates
306  Vector<double> x(2);
307  for(unsigned i=0;i<2;i++)
308  {
309  x[i]=this->interpolated_x(s,i);
310  outfile << x[i] << " ";
311  }
312 
313  // Veloc
314  outfile << veloc[0] << " "
315  << veloc[1] << " ";
316 
317  // Smoothed vorticity
318  outfile << vort << " ";
319 
320  // Smoothed vorticity derivatives (d/dx, d/dy, d^2/dx^2, d^2/dxdy, d^2/dy^2
321  // d^3/dx^3, d^3/dx^2dy, d^3/dxdy^2, d^3/dy^3
322  // du/dx, du/dy, dv/dx, dv/dy
323  outfile << dvort_dx[0] << " "
324  << dvort_dx[1] << " "
325  << dvort_dxdy[0] << " "
326  << dvort_dxdy[1] << " "
327  << dvort_dxdy[2] << " "
328  << dvort_dxdxdy[0] << " "
329  << dvort_dxdxdy[1] << " "
330  << dvort_dxdxdy[2] << " "
331  << dvort_dxdxdy[3] << " "
332  << dveloc_dx[0] << " "
333  << dveloc_dx[1] << " "
334  << dveloc_dx[2] << " "
335  << dveloc_dx[3] << " ";
336 
337  outfile << std::endl;
338  }
339 
340  // Write tecplot footer (e.g. FE connectivity lists)
341  this->write_tecplot_zone_footer(outfile,nplot);
342  }
343 
344 
345 
346 
347 
348 
349 
350  /// \short Overloaded output fct: Output veloc, pressure,
351  /// smoothed vorticity
352  void output(std::ostream &outfile, const unsigned &nplot)
353  {
354  //Vector of local coordinates
355  Vector<double> s(2);
356 
357  // Shape functions
358  unsigned n_node = this->nnode();
359  Shape psif(n_node);
360 
361  // Tecplot header info
362  outfile << this->tecplot_zone_string(nplot);
363 
364  // Loop over plot points
365  unsigned num_plot_points=this->nplot_points(nplot);
366  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
367  {
368  // Get local coordinates of plot point
369  this->get_s_plot(iplot,nplot,s);
370 
371  // Get shape fct
372  this->shape(s,psif);
373 
374  // Coordinates
375  Vector<double> x(2);
376  for(unsigned i=0;i<2;i++)
377  {
378  x[i]=this->interpolated_x(s,i);
379  outfile << x[i] << " ";
380  }
381 
382  // Veloc
383  Vector<double> veloc(2);
384  this->interpolated_u_nst(s,veloc);
385  outfile << veloc[0] << " "
386  << veloc[1] << " ";
387 
388  // Pressure
389  outfile << this->interpolated_p_nst(s) << " ";
390 
391  // Smoothed vorticity
392  double smoothed_vort=0.0;
393  for(unsigned l=0;l<n_node;l++)
394  {
395  smoothed_vort += this->nodal_value(l,Smoothed_vorticity_index)*psif[l];
396  }
397  outfile << smoothed_vort << " ";
398 
399  // Smoothed vorticity derivatives (d/dx, d/dy, d^2/dx^2, d^2/dxdy, d^2/dy^2
400  // d^3/dx^3, d^3/dx^2dy, d^3/dxdy^2, d^3/dy^3,
401  // du/dx, du/dy, dv/dx, dv/dy
402  for (unsigned i=1;i<14;i++) // hierher was 10
403  {
404  double smoothed_vort_deriv=0.0;
405  for(unsigned l=0;l<n_node;l++)
406  {
407  smoothed_vort_deriv+=
408  this->nodal_value(l,Smoothed_vorticity_index+i)*psif[l];
409  }
410  outfile << smoothed_vort_deriv << " ";
411  }
412 
413  outfile << std::endl;
414  }
415 
416  // Write tecplot footer (e.g. FE connectivity lists)
417  this->write_tecplot_zone_footer(outfile,nplot);
418  }
419 
420 
421  /// Get raw derivative of velocity
423  Vector<double>& dveloc_dx) const
424  {
425  //Find out how many nodes there are
426  unsigned n_node = this->nnode();
427 
428  //Set up memory for the shape functions
429  Shape psif(n_node);
430  DShape dpsifdx(n_node,2);
431 
432  //Call the derivatives of the shape and test functions
433  this->dshape_eulerian(s,psif,dpsifdx);
434 
435  //Initialise to zero
436  for(unsigned j=0;j<4;j++)
437  {
438  dveloc_dx[j] = 0.0;
439  }
440 
441  // Loop over nodes
442  for(unsigned l=0;l<n_node;l++)
443  {
444  //Loop over derivative directions
445  for(unsigned j=0;j<2;j++)
446  {
447  dveloc_dx[j] += this->nodal_value(l,0)*dpsifdx(l,j);
448  dveloc_dx[j+2] += this->nodal_value(l,1)*dpsifdx(l,j);
449  }
450  }
451  }
452 
453  /// Get raw derivative of smoothed vorticity
455  Vector<double>& dvorticity_dx) const
456  {
457  //Find out how many nodes there are
458  unsigned n_node = this->nnode();
459 
460  //Set up memory for the shape functions
461  Shape psif(n_node);
462  DShape dpsifdx(n_node,2);
463 
464  //Call the derivatives of the shape and test functions
465  this->dshape_eulerian(s,psif,dpsifdx);
466 
467  //Initialise to zero
468  for(unsigned j=0;j<2;j++)
469  {
470  dvorticity_dx[j] = 0.0;
471  }
472 
473  // Loop over nodes
474  for(unsigned l=0;l<n_node;l++)
475  {
476  //Loop over derivative directions
477  for(unsigned j=0;j<2;j++)
478  {
479  dvorticity_dx[j] += this->nodal_value(l,Smoothed_vorticity_index)*
480  dpsifdx(l,j);
481  }
482  }
483  }
484 
485  /// Get raw derivative of smoothed derivative vorticity
487  Vector<double>& dvorticity_dxdy) const
488  {
489  //Find out how many nodes there are
490  unsigned n_node = this->nnode();
491 
492  //Set up memory for the shape functions
493  Shape psif(n_node);
494  DShape dpsifdx(n_node,2);
495 
496  //Call the derivatives of the shape and test functions
497  this->dshape_eulerian(s,psif,dpsifdx);
498 
499  //Initialise to zero
500  for(unsigned j=0;j<3;j++)
501  {
502  dvorticity_dxdy[j] = 0.0;
503  }
504 
505  // Loop over nodes
506  for(unsigned l=0;l<n_node;l++)
507  {
508  //Loop over derivative directions to obtain xx and xy derivatives
509  for(unsigned j=0;j<2;j++)
510  {
511  dvorticity_dxdy[j] += this->nodal_value(l,Smoothed_vorticity_index+1)*
512  dpsifdx(l,j);
513  }
514  //Calcutation of yy derivative
515  dvorticity_dxdy[2] += this->nodal_value(l,Smoothed_vorticity_index+2)*
516  dpsifdx(l,1);
517  }
518  }
519 
520 
521  /// Get raw derivative of smoothed derivative vorticity
522  /// [0]: d^3/dx^3, [0]: d^3/dx^2dy, [0]: d^3/dxdy^2, [0]: d^3/dy^3,
524  Vector<double>& dvorticity_dxdy) const
525  {
526 
527  //Find out how many nodes there are
528  unsigned n_node = this->nnode();
529 
530  //Set up memory for the shape functions
531  Shape psif(n_node);
532  DShape dpsifdx(n_node,2);
533 
534  //Call the derivatives of the shape and test functions
535  this->dshape_eulerian(s,psif,dpsifdx);
536 
537  //Initialise to zero
538  for(unsigned j=0;j<4;j++)
539  {
540  dvorticity_dxdy[j] = 0.0;
541  }
542 
543  // Loop over nodes
544  for(unsigned l=0;l<n_node;l++)
545  {
546  // d^3/dx^3 = d/dx \overline{d^2/dx^2}
547  dvorticity_dxdy[0] += this->nodal_value(l,Smoothed_vorticity_index+3)*
548  dpsifdx(l,0);
549 
550  // d^3/dx^2dy = d/dx \overline{d^2/dxdy}
551  dvorticity_dxdy[1] += this->nodal_value(l,Smoothed_vorticity_index+4)*
552  dpsifdx(l,0);
553 
554  // d^3/dxdy^2 = d/dy \overline{d^2/dxdy}
555  dvorticity_dxdy[2] += this->nodal_value(l,Smoothed_vorticity_index+4)*
556  dpsifdx(l,1);
557 
558  // d^3/dy^3 = d/dy \overline{d^2/dy^2}
559  dvorticity_dxdy[3] += this->nodal_value(l,Smoothed_vorticity_index+5)*
560  dpsifdx(l,1);
561 
562  }
563  }
564 
565 
566  /// \short Compute the element's contribution to the (squared) L2 norm
567  /// of the difference between exact and smoothed vorticity. i=0: do
568  /// vorticity itself; i>0: derivs
569  double vorticity_error_squared(const unsigned& i)
570  {
571  double norm_squared=0.0;
572 
573  //Find out how many nodes there are
574  unsigned n_node = this->nnode();
575 
576  //Set up memory for the shape functions
577  Shape psif(n_node);
578  DShape dpsifdx(n_node,2);
579 
580  //Number of integration points
581  unsigned n_intpt = this->integral_pt()->nweight();
582 
583  //Set the Vector to hold local coordinates
584  Vector<double> s(2);
585  Vector<double> x(2);
586 
587  //Loop over the integration points
588  for(unsigned ipt=0;ipt<n_intpt;ipt++)
589  {
590  //Assign values of s
591  for(unsigned ii=0;ii<2;ii++)
592  {
593  x[ii]=0.0;
594  s[ii] = this->integral_pt()->knot(ipt,ii);
595  }
596 
597  //Get the integral weight
598  double w = this->integral_pt()->weight(ipt);
599 
600  //Call the derivatives of the shape and test functions
601  double J = this->dshape_eulerian(s,psif,dpsifdx);
602 
603  //Premultiply the weights and the Jacobian
604  double W = w*J;
605 
606  // Smoothed vorticity
607  double smoothed_vort=0.0;
608  for(unsigned l=0;l<n_node;l++)
609  {
610  Node* nod_pt=this->node_pt(l);
611  smoothed_vort += this->nodal_value(l,Smoothed_vorticity_index+i)*psif[l];
612  for (unsigned ii=0;ii<2;ii++)
613  {
614  x[ii]+=nod_pt->x(ii)*psif[l];
615  }
616  }
617 
618  // Synthetic quantities
619  double synth_vort=0.0;
620  Vector<double> synth_dvort_dx(2,0.0);
621  Vector<double> synth_dvort_dxdy(3,0.0);
622  Vector<double> synth_dvort_dxdxdy(4,0.0);
623  Vector<double> synth_dveloc_dx(4,0.0);
624 
625  if (Exact_vorticity_fct_pt!=0)
626  {
627  Exact_vorticity_fct_pt(x,synth_vort,synth_dvort_dx,
628  synth_dvort_dxdy,synth_dvort_dxdxdy,
629  synth_dveloc_dx);
630  }
631  double synth_quantity=synth_vort;
632  if ((i==1)||(i==2))
633  {
634  synth_quantity=synth_dvort_dx[i-1];
635  }
636  if ((i==3)||(i==4)||(i==5))
637  {
638  synth_quantity=synth_dvort_dxdy[i-3];
639  }
640  if ((i==6)||(i==7)||(i==8)||(i==9))
641  {
642  synth_quantity=synth_dvort_dxdxdy[i-6];
643  }
644  if ((i==10)||(i==11)||(i==12)||(i==13))
645  {
646  synth_quantity=synth_dveloc_dx[i-10];
647  }
648  // Add squared difference
649  norm_squared+=pow(smoothed_vort-synth_quantity,2)*W;
650  }
651 
652  return norm_squared;
653  }
654 
655 
656  /// \short Compute smoothed vorticity and its derivatives
658  Vector<double>& veloc,
659  double& vort,
660  Vector<double>& dvort_dx,
661  Vector<double>& dvort_dxdy,
662  Vector<double>& dvort_dxdxdy,
663  Vector<double>& dveloc_dx)
664  {
665  // Shape functions
666  unsigned n_node = this->nnode();
667  Shape psif(n_node);
668  this->shape(s,psif);
669 
670  // Smoothed vorticity
671  vort=0.0;
672  veloc[0]=0.0;
673  veloc[1]=0.0;
674  for(unsigned l=0;l<n_node;l++)
675  {
676  vort += this->nodal_value(l,Smoothed_vorticity_index)*psif[l];
677  veloc[0]+=this->nodal_value(l,0)*psif[l];
678  veloc[1]+=this->nodal_value(l,1)*psif[l];
679  }
680 
681  // Smoothed vorticity derivatives (d/dx, d/dy, d^2/dx^2, d^2/dxdy, d^2/dy^2,
682  // d^3/dx^3, d^3/dx^2dy, d^3/dxdy^2,d^3,dy^3
683  // du/dx, du/dy, dv/dx, dv/dy
684  for (unsigned i=1;i<14;i++) // hierher used to be 10
685  {
686  double smoothed_vort_deriv=0.0;
687  for(unsigned l=0;l<n_node;l++)
688  {
689  smoothed_vort_deriv+=
690  this->nodal_value(l,Smoothed_vorticity_index+i)*psif[l];
691  }
692  switch (i)
693  {
694  case 1:
695  dvort_dx[0]=smoothed_vort_deriv;
696  break;
697  case 2:
698  dvort_dx[1]=smoothed_vort_deriv;
699  break;
700  case 3:
701  dvort_dxdy[0]=smoothed_vort_deriv;
702  break;
703  case 4:
704  dvort_dxdy[1]=smoothed_vort_deriv;
705  break;
706  case 5:
707  dvort_dxdy[2]=smoothed_vort_deriv;
708  break;
709  case 6:
710  dvort_dxdxdy[0]=smoothed_vort_deriv;
711  break;
712  case 7:
713  dvort_dxdxdy[1]=smoothed_vort_deriv;
714  break;
715  case 8:
716  dvort_dxdxdy[2]=smoothed_vort_deriv;
717  break;
718  case 9:
719  dvort_dxdxdy[3]=smoothed_vort_deriv;
720  break;
721  case 10:
722  dveloc_dx[0]=smoothed_vort_deriv;
723  break;
724  case 11:
725  dveloc_dx[1]=smoothed_vort_deriv;
726  break;
727  case 12:
728  dveloc_dx[2]=smoothed_vort_deriv;
729  break;
730  case 13:
731  dveloc_dx[3]=smoothed_vort_deriv;
732  break;
733  default:
734  oomph_info << "never get here\n";
735  abort();
736  break;
737  }
738  }
739  }
740 
741  private:
742 
743  /// Index of smoothed vorticity -- followed by derivatives
745 
746  /// \short Number of values per field: u,v,p,omega,d/dx,d/dy,
747  /// d^2/dx^2,d^2/dxdy, d^2/dy^2,
748  // du/dx, du/dy, dv/dx, dv/dy
750 
751  /// Pointer to fct that specifies exact vorticity and
752  /// derivs (for validation).
754 
755 
756 };
757 
758 
759 } // end namespace extension
760 
761 
762 
763 
764 
765 
766 //////////////////////////////////////////////////////////////////////
767 //////////////////////////////////////////////////////////////////////
768 //////////////////////////////////////////////////////////////////////
769 
770 
771 
772 
773 //========================================================
774 /// Smoother for vorticity in 2D
775 //========================================================
776 template<class ELEMENT>
778 {
779  public:
780 
781 
782  /// Constructor: Set order of recovery shape functions
783  VorticitySmoother(const unsigned& recovery_order) :
784  Recovery_order(recovery_order)
785  {}
786 
787  /// Broken copy constructor
789  {
790  BrokenCopy::broken_copy("VorticitySmoother");
791  }
792 
793  /// Broken assignment operator
795  {
796  BrokenCopy::broken_assign("VorticitySmoother");
797  }
798 
799  /// Empty virtual destructor
800  virtual ~VorticitySmoother(){}
801 
802  /// Access function for order of recovery polynomials
803  unsigned& recovery_order() {return Recovery_order;}
804 
805  /// Recovery shape functions as functions of the global, Eulerian
806  /// coordinate x of dimension dim.
807  /// The recovery shape functions are complete polynomials of
808  /// the order specified by Recovery_order.
809  void shape_rec(const Vector<double>& x,
810  Vector<double>& psi_r)
811  {
812  std::ostringstream error_stream;
813 
814  /// Find order of recovery shape functions
815  switch(recovery_order())
816  {
817  case 1:
818 
819  // Complete linear polynomial in 2D:
820  psi_r[0]=1.0;
821  psi_r[1]=x[0];
822  psi_r[2]=x[1];
823  break;
824 
825  case 2:
826 
827  // Complete quadratic polynomial in 2D:
828  psi_r[0]=1.0;
829  psi_r[1]=x[0];
830  psi_r[2]=x[1];
831  psi_r[3]=x[0]*x[0];
832  psi_r[4]=x[0]*x[1];
833  psi_r[5]=x[1]*x[1];
834  break;
835 
836  case 3:
837 
838  // Complete cubic polynomial in 2D:
839  psi_r[0]=1.0;
840  psi_r[1]=x[0];
841  psi_r[2]=x[1];
842  psi_r[3]=x[0]*x[0];
843  psi_r[4]=x[0]*x[1];
844  psi_r[5]=x[1]*x[1];
845  psi_r[6]=x[0]*x[0]*x[0];
846  psi_r[7]=x[0]*x[0]*x[1];
847  psi_r[8]=x[0]*x[1]*x[1];
848  psi_r[9]=x[1]*x[1]*x[1];
849  break;
850 
851  default:
852 
853  error_stream
854  << "Recovery shape functions for recovery order "
855  << recovery_order() << " haven't yet been implemented for 2D"
856  << std::endl;
857 
858  throw OomphLibError(error_stream.str(),
859  OOMPH_CURRENT_FUNCTION,
860  OOMPH_EXCEPTION_LOCATION);
861  }
862  }
863 
864 
865 
866 /// Integation scheme associated with the recovery shape functions
867 /// must be of sufficiently high order to integrate the mass matrix
868 /// associated with the recovery shape functions. The argument
869 /// is the dimension of the elements.
870 /// The integration is performed locally over the elements, so the
871 /// integration scheme does depend on the geometry of the element.
872 /// The type of element is specified by the boolean which is
873 /// true if elements in the patch are QElements and false if they are
874 /// TElements (will need change if we ever have other element types)
875  Integral* integral_rec(const bool &is_q_mesh)
876  {
877  std::ostringstream error_stream;
878 
879  // 2D:
880 
881  /// Find order of recovery shape functions
882  switch(recovery_order())
883  {
884  case 1:
885 
886  // Complete linear polynomial in 2D:
887  if(is_q_mesh) {return(new Gauss<2,2>);}
888  else {return(new TGauss<2,2>);}
889  break;
890 
891  case 2:
892 
893  // Complete quadratic polynomial in 2D:
894  if(is_q_mesh) {return(new Gauss<2,3>);}
895  else {return(new TGauss<2,3>);}
896  break;
897 
898  case 3:
899 
900  // Complete cubic polynomial in 2D:
901  if(is_q_mesh) {return(new Gauss<2,4>);}
902  else {return(new TGauss<2,4>);}
903  break;
904 
905  default:
906 
907  error_stream
908  << "Recovery shape functions for recovery order "
909  << recovery_order() << " haven't yet been implemented for 2D"
910  << std::endl;
911 
912  throw OomphLibError(error_stream.str(),
913  OOMPH_CURRENT_FUNCTION,
914  OOMPH_EXCEPTION_LOCATION);
915  }
916 
917 
918  //Dummy return (never get here)
919  return 0;
920  }
921 
922 
923  /// \short Setup patches: For each vertex node pointed to by nod_pt,
924  /// adjacent_elements_pt[nod_pt] contains the pointer to the vector that
925  /// contains the pointers to the elements that the node is part of.
926  /// Also returns a Vector of vertex nodes for use in get_element_errors.
927  void setup_patches(Mesh*& mesh_pt,
928  std::map<Node*,Vector<ELEMENT*>*>& adjacent_elements_pt,
929  Vector<Node*>& vertex_node_pt)
930  {
931 
932  // Clear: hierher should we do this in Z2 as well?
933  adjacent_elements_pt.clear();
934 
935  // Auxiliary map that contains element-adjacency for ALL nodes
936  std::map<Node*,Vector<ELEMENT*>*> aux_adjacent_elements_pt;
937 
938 #ifdef PARANOID
939  // Check if all elements request the same recovery order
940  unsigned ndisagree=0;
941 #endif
942 
943  // Loop over all elements to setup adjacency for all nodes.
944  // Need to do this because midside nodes can be corner nodes for
945  // adjacent smaller elements! Admittedly, the inclusion of interior
946  // nodes is wasteful...
947  unsigned nelem=mesh_pt->nelement();
948  for (unsigned e=0;e<nelem;e++)
949  {
950  ELEMENT* el_pt=
951  dynamic_cast<ELEMENT*>(mesh_pt->element_pt(e));
952 
953 #ifdef PARANOID
954  // Check if all elements request the same recovery order
955  if (el_pt->nrecovery_order()!=Recovery_order){ndisagree++;}
956 #endif
957 
958  // Loop all nodes in element
959  unsigned nnod=el_pt->nnode();
960  for (unsigned n=0;n<nnod;n++)
961  {
962  Node* nod_pt=el_pt->node_pt(n);
963 
964  // Has this node been considered before?
965  if (aux_adjacent_elements_pt[nod_pt]==0)
966  {
967  // Create Vector of pointers to its adjacent elements
968  aux_adjacent_elements_pt[nod_pt]=new
969  Vector<ELEMENT*>;
970  }
971 
972  // Add pointer to adjacent element
973  (*aux_adjacent_elements_pt[nod_pt]).push_back(el_pt);
974  }
975  } // end element loop
976 
977 #ifdef PARANOID
978  // Check if all elements request the same recovery order
979  if (ndisagree!=0)
980  {
981  oomph_info
982  << "\n\n========================================================\n";
983  oomph_info << "WARNING: " << std::endl;
984  oomph_info << ndisagree << " out of " << mesh_pt->nelement()
985  << " elements\n";
986  oomph_info << "have different preferences for the order of the recovery\n";
987  oomph_info << "shape functions. We are using: Recovery_order="
988  << Recovery_order << std::endl;
989  oomph_info
990  << "========================================================\n\n";
991  }
992 #endif
993 
994  //Loop over all elements, extract adjacency for corner nodes only
995  nelem=mesh_pt->nelement();
996  for (unsigned e=0;e<nelem;e++)
997  {
998  ELEMENT* el_pt=
999  dynamic_cast<ELEMENT*>(mesh_pt->element_pt(e));
1000 
1001  // Loop over corner nodes
1002  unsigned n_node=el_pt->nvertex_node();
1003  for (unsigned n=0;n<n_node;n++)
1004  {
1005  Node* nod_pt=el_pt->vertex_node_pt(n);
1006 
1007  // Has this node been considered before?
1008  if (adjacent_elements_pt[nod_pt]==0)
1009  {
1010 
1011  // Add the node pointer to the vertex node container
1012  vertex_node_pt.push_back(nod_pt);
1013 
1014  // Create Vector of pointers to its adjacent elements
1015  adjacent_elements_pt[nod_pt]=new
1016  Vector<ELEMENT*>;
1017 
1018  // Copy across:
1019  unsigned nel=(*aux_adjacent_elements_pt[nod_pt]).size();
1020  for (unsigned e=0;e<nel;e++)
1021  {
1022  (*adjacent_elements_pt[nod_pt]).push_back(
1023  (*aux_adjacent_elements_pt[nod_pt])[e]);
1024  }
1025  }
1026  }
1027 
1028  } // end of loop over elements
1029 
1030 
1031  // Cleanup
1032  for (typename std::map<Node*,Vector<ELEMENT*>*>::iterator it=
1033  aux_adjacent_elements_pt.begin();
1034  it!=aux_adjacent_elements_pt.end();it++)
1035  {
1036  delete it->second;
1037  }
1038 
1039  }
1040 
1041 
1042 
1043  /// Given the vector of elements that make up a patch, compute
1044  /// the vector of recovered vorticity coefficients and return
1045  /// a pointer to it. n_deriv indicates which derivative of the
1046  /// vorticity is supposed to be smoothed: 0: zeroth (i.e. vorticity
1047  /// itself; 1: d/dx; 2: d/dy; 3: d^2/dx^2; 4: d^2/dxdy 5: d^2/dy^2
1048  /// 6: d^3/dx^3, 7: d^3/dx^2dy, 8: d^3/dxdy^2, 9: d^3/dy^3,
1049  /// 10: du/dx, 11: du/dy, 12: dv/dx, 13: dv/dy
1051  const Vector<ELEMENT*>& patch_el_pt,
1052  const unsigned& num_recovery_terms,
1053  Vector<double>*& recovered_vorticity_coefficient_pt,
1054  unsigned& n_deriv)
1055  {
1056 
1057  // Create/initialise matrix for linear system
1058  DenseDoubleMatrix recovery_mat(num_recovery_terms,num_recovery_terms,0.0);
1059 
1060  // Ceate/initialise vector for RHS
1061  Vector<double> rhs(num_recovery_terms,0.0);
1062 
1063  //Create a new integration scheme based on the recovery order
1064  //in the elements
1065  //Need to find the type of the element, default is to assume a quad
1066  bool is_q_mesh=true;
1067 
1068  //If we can dynamic cast to the TElementBase, then it's a triangle/tet
1069  //Note that I'm assuming that all elements are of the same geometry, but
1070  //if they weren't we could adapt...
1071  if(dynamic_cast<TElementBase*>(patch_el_pt[0])) {is_q_mesh=false;}
1072 
1073  Integral* const integ_pt = this->integral_rec(is_q_mesh);
1074 
1075  //Loop over all elements in patch to assemble linear system
1076  unsigned nelem=patch_el_pt.size();
1077  for (unsigned e=0;e<nelem;e++)
1078  {
1079  // Get pointer to element
1080  ELEMENT* const el_pt=patch_el_pt[e];
1081 
1082  // Create storage for the recovery shape function values
1083  Vector<double> psi_r(num_recovery_terms);
1084 
1085  //Create vector to hold local coordinates
1086  Vector<double> s(2);
1087 
1088  //Loop over the integration points
1089  unsigned n_intpt = integ_pt->nweight();
1090  for(unsigned ipt=0;ipt<n_intpt;ipt++)
1091  {
1092  //Assign values of s, the local coordinate
1093  for(unsigned i=0;i<2;i++)
1094  {
1095  s[i] = integ_pt->knot(ipt,i);
1096  }
1097 
1098  //Get the integral weight
1099  double w = integ_pt->weight(ipt);
1100 
1101  //Jacobian of mapping
1102  double J = el_pt->J_eulerian(s);
1103 
1104  // Interpolate the global (Eulerian) coordinate
1105  Vector<double> x(2);
1106  el_pt->interpolated_x(s,x);
1107 
1108  // Premultiply the weights and the Jacobian
1109  // and the geometric jacobian weight (used in axisymmetric
1110  // and spherical coordinate systems) -- hierher really fct of x?
1111  // probably yes, actually).
1112  double W = w*J*(el_pt->geometric_jacobian(x));
1113 
1114  // Recovery shape functions at global (Eulerian) coordinate
1115  shape_rec(x,psi_r);
1116 
1117  // Get FE estimates for vorticity:
1118  Vector<double> vorticity(1);
1119  if(n_deriv==0)
1120  {
1121  el_pt->get_vorticity(s,vorticity);
1122  }
1123 
1124  // Get FE estimates for deriv of vorticity:
1125  Vector<double> deriv_vorticity(2);
1126  Vector<double> second_deriv_vorticity(3);
1127  Vector<double> third_deriv_vorticity(4);
1128  Vector<double> deriv_velocity(4);
1129  if((n_deriv==1)||(n_deriv==2))
1130  {
1131  el_pt->get_raw_vorticity_deriv(s,deriv_vorticity);
1132  }
1133  // Get FE estimates for second deriv of vorticity:
1134  else if((n_deriv==3)||(n_deriv==4)||(n_deriv==5))
1135  {
1136  el_pt->get_raw_vorticity_second_deriv(s,second_deriv_vorticity);
1137  }
1138  // Get FE estimates for third deriv of vorticity:
1139  else if((n_deriv==6)||(n_deriv==7)||(n_deriv==8)||(n_deriv==9))
1140  {
1141  el_pt->get_raw_vorticity_third_deriv(s,third_deriv_vorticity);
1142  }
1143  // Get FE estimates for derivs of velocity
1144  else if((n_deriv==10)||(n_deriv==11)||(n_deriv==12)||(n_deriv==13))
1145  {
1146  el_pt->get_raw_velocity_deriv(s,deriv_velocity);
1147  }
1148 
1149  // Add elemental RHSs and recovery matrix to global versions
1150  //----------------------------------------------------------
1151  if(n_deriv==0)
1152  {
1153  // RHS
1154  // Loop over the nodes for the test functions
1155  for(unsigned l=0;l<num_recovery_terms;l++)
1156  {
1157  rhs[l] += vorticity[0]*psi_r[l]*W;
1158  }
1159  }
1160  // Get x-derivative of vorticity
1161  else if(n_deriv==1)
1162  {
1163  // RHS
1164  // Loop over the nodes for the test functions
1165  for(unsigned l=0;l<num_recovery_terms;l++)
1166  {
1167  rhs[l] += deriv_vorticity[0]*psi_r[l]*W;
1168  }
1169  }
1170  // Get y-derivative of vorticity
1171  else if(n_deriv==2)
1172  {
1173  // RHS
1174  // Loop over the nodes for the test functions
1175  for(unsigned l=0;l<num_recovery_terms;l++)
1176  {
1177  rhs[l] += deriv_vorticity[1]*psi_r[l]*W;
1178  }
1179  }
1180  // Get xx-derivative of vorticity
1181  else if(n_deriv==3)
1182  {
1183  // RHS
1184  // Loop over the nodes for the test functions
1185  for(unsigned l=0;l<num_recovery_terms;l++)
1186  {
1187  rhs[l] += second_deriv_vorticity[0]*psi_r[l]*W;
1188  }
1189  }
1190  // Get xy-derivative of vorticity
1191  else if(n_deriv==4)
1192  {
1193  // RHS
1194  // Loop over the nodes for the test functions
1195  for(unsigned l=0;l<num_recovery_terms;l++)
1196  {
1197  rhs[l] += second_deriv_vorticity[1]*psi_r[l]*W;
1198  }
1199  }
1200  // Get yy-derivative of vorticity
1201  else if(n_deriv==5)
1202  {
1203  // RHS
1204  // Loop over the nodes for the test functions
1205  for(unsigned l=0;l<num_recovery_terms;l++)
1206  {
1207  rhs[l] += second_deriv_vorticity[2]*psi_r[l]*W;
1208  }
1209  }
1210  // Get xxx-derivative of vorticity
1211  else if(n_deriv==6)
1212  {
1213  // RHS
1214  // Loop over the nodes for the test functions
1215  for(unsigned l=0;l<num_recovery_terms;l++)
1216  {
1217  rhs[l] += third_deriv_vorticity[0]*psi_r[l]*W;
1218  }
1219  }
1220  // Get xxy-derivative of vorticity
1221  else if(n_deriv==7)
1222  {
1223  // RHS
1224  // Loop over the nodes for the test functions
1225  for(unsigned l=0;l<num_recovery_terms;l++)
1226  {
1227  rhs[l] += third_deriv_vorticity[1]*psi_r[l]*W;
1228  }
1229  }
1230  // Get xxy-derivative of vorticity
1231  else if(n_deriv==8)
1232  {
1233  // RHS
1234  // Loop over the nodes for the test functions
1235  for(unsigned l=0;l<num_recovery_terms;l++)
1236  {
1237  rhs[l] += third_deriv_vorticity[2]*psi_r[l]*W;
1238  }
1239  }
1240  // Get xxy-derivative of vorticity
1241  else if(n_deriv==9)
1242  {
1243  // RHS
1244  // Loop over the nodes for the test functions
1245  for(unsigned l=0;l<num_recovery_terms;l++)
1246  {
1247  rhs[l] += third_deriv_vorticity[3]*psi_r[l]*W;
1248  }
1249  }
1250  // Get x-derivative of u-velocity
1251  else if(n_deriv==10)
1252  {
1253  // RHS
1254  // Loop over the nodes for the test functions
1255  for(unsigned l=0;l<num_recovery_terms;l++)
1256  {
1257  rhs[l] += deriv_velocity[0]*psi_r[l]*W;
1258  }
1259  }
1260  // Get y-derivative of u-velocity
1261  else if(n_deriv==11)
1262  {
1263  // RHS
1264  // Loop over the nodes for the test functions
1265  for(unsigned l=0;l<num_recovery_terms;l++)
1266  {
1267  rhs[l] += deriv_velocity[1]*psi_r[l]*W;
1268  }
1269  }
1270  // Get x-derivative of v-velocity
1271  else if(n_deriv==12)
1272  {
1273  // RHS
1274  // Loop over the nodes for the test functions
1275  for(unsigned l=0;l<num_recovery_terms;l++)
1276  {
1277  rhs[l] += deriv_velocity[2]*psi_r[l]*W;
1278  }
1279  }
1280  // Get y-derivative of v-velocity
1281  else if(n_deriv==13)
1282  {
1283  // RHS
1284  // Loop over the nodes for the test functions
1285  for(unsigned l=0;l<num_recovery_terms;l++)
1286  {
1287  rhs[l] += deriv_velocity[3]*psi_r[l]*W;
1288  }
1289  }
1290  else
1291  {
1292  oomph_info << "Never get here\n";
1293  abort();
1294  }
1295 
1296 
1297  // Loop over the nodes for the test functions
1298  for(unsigned l=0;l<num_recovery_terms;l++)
1299  {
1300  //Loop over the nodes for the variables
1301  for(unsigned l2=0;l2<num_recovery_terms;l2++)
1302  {
1303  //Add contribution to recovery matrix
1304  recovery_mat(l,l2)+=psi_r[l]*psi_r[l2]*W;
1305  }
1306  }
1307  }
1308  } // End of loop over elements that make up patch.
1309 
1310  //Delete the integration scheme
1311  delete integ_pt;
1312 
1313  // Linear system is now assembled: Solve recovery system
1314 
1315  // LU decompose the recovery matrix
1316  recovery_mat.ludecompose();
1317 
1318  // Back-substitute
1319  recovery_mat.lubksub(rhs);
1320 
1321  // Now create a matrix to store the vorticity recovery coefficients.
1322  // Pointer to this matrix will be returned.
1323  recovered_vorticity_coefficient_pt =
1324  new Vector<double>(num_recovery_terms);
1325 
1326  // Copy coefficients
1327  for (unsigned icoeff=0;icoeff<num_recovery_terms;icoeff++)
1328  {
1329  (*recovered_vorticity_coefficient_pt)[icoeff]=rhs[icoeff];
1330  }
1331 
1332 }
1333 
1334  // Get the recovery order
1335  unsigned nrecovery_order() const
1336  {
1337  switch(Recovery_order)
1338  {
1339  case 1:
1340 
1341  // Linear recovery shape functions
1342  //--------------------------------
1343  return 3; // 1, x, y
1344  break;
1345 
1346 
1347  case 2:
1348 
1349  // Quadratic recovery shape functions
1350  //-----------------------------------
1351  return 6; // 1, x, y, x^2, xy, y^2
1352  break;
1353 
1354  case 3:
1355 
1356  // Cubic recovery shape functions
1357  //--------------------------------
1358  return 10; // 1, x, y, x^2, xy, y^2, x^3, y^3, x^2 y, x y^2
1359  break;
1360 
1361  default:
1362 
1363  // Any other recovery order?
1364  //--------------------------
1365  std::ostringstream error_stream;
1366  error_stream
1367  << "Wrong Recovery_order " << Recovery_order << std::endl;
1368 
1369  throw OomphLibError(error_stream.str(),
1370  OOMPH_CURRENT_FUNCTION,
1371  OOMPH_EXCEPTION_LOCATION);
1372  }
1373  }
1374 
1375 
1376  /// Recover vorticity from patches
1377  void recover_vorticity(Mesh* mesh_pt)
1378  {
1379  DocInfo doc_info;
1380  doc_info.disable_doc();
1381  recover_vorticity(mesh_pt,doc_info);
1382  }
1383 
1384 
1385 
1386  /// Recover vorticity from patches -- output intermediate steps
1387  /// to directory specified by DocInfo object
1388  void recover_vorticity(Mesh* mesh_pt, DocInfo& doc_info)
1389  {
1390 
1391  double t_start=TimingHelpers::timer();
1392 
1393  // Make patches
1394  //-------------
1395  std::map<Node*,Vector<ELEMENT*>*> adjacent_elements_pt;
1396  Vector<Node*> vertex_node_pt;
1397  setup_patches(mesh_pt,
1398  adjacent_elements_pt,
1399  vertex_node_pt);
1400 
1401  // Determine number of coefficients for expansion of recovered vorticity
1402  // Use complete polynomial of given order for recovery
1403  unsigned num_recovery_terms=nrecovery_order();
1404 
1405  // hierher get from element
1406  unsigned Smoothed_vorticity_index=3;
1407 
1408  // Counter for averaging of recovered vorticity and its derivatives
1409  map<Node* , unsigned> count;
1410 
1411  // Loop over derivatives
1412  for (unsigned deriv=0;deriv<14;deriv++) // hierher used to be 10
1413  {
1414 
1415  // Storage for accumulated nodal vorticity (used to compute
1416  // nodal averages)
1417  map<Node*, double> averaged_recovered_vort;
1418 
1419  //Calculation of vorticity
1420  //------------------------
1421 
1422  // Do patch recovery
1423  // unsigned counter=0;
1424  for (typename std::map<Node*,Vector<ELEMENT*>*>::iterator it=
1425  adjacent_elements_pt.begin();
1426  it!=adjacent_elements_pt.end();it++)
1427  {
1428 
1429  // Setup smoothed vorticity field for patches
1430  Vector<double>* recovered_vorticity_coefficient_pt;
1431  get_recovered_vorticity_in_patch(*(it->second),
1432  num_recovery_terms,
1433  recovered_vorticity_coefficient_pt,
1434  deriv);
1435 
1436  // Now get the nodal average of the recovered vorticity
1437  // (nodes are generally part of multiple patches)
1438 
1439  //Loop over all elements to get recovered vorticity
1440  unsigned nelem=(*(it->second)).size();
1441  for (unsigned e=0;e<nelem;e++)
1442  {
1443  // Get pointer to element
1444  ELEMENT* const el_pt=(*(it->second))[e];
1445 
1446  // Get the number of nodes by element
1447  unsigned nnode_el=el_pt->nnode();
1448  for(unsigned j=0;j<nnode_el;j++)
1449  {
1450  //Get local coordinates
1451  Vector<double> s(2);
1452  Node* nod_pt=el_pt->node_pt(j);
1453  el_pt->local_coordinate_of_node(j,s);
1454 
1455  // Interpolate the global (Eulerian) coordinate
1456  Vector<double> x(2);
1457  el_pt->interpolated_x(s,x);
1458 
1459  // Recovery shape functions at global (Eulerian) coordinate
1460  Vector<double> psi_r(num_recovery_terms);
1461  shape_rec(x,psi_r);
1462 
1463  // Assemble recovered vorticity
1464  double recovered_vort=0.0;
1465  for (unsigned i=0;i<num_recovery_terms;i++)
1466  {
1467  recovered_vort+=(*recovered_vorticity_coefficient_pt)[i]*psi_r[i];
1468  }
1469 
1470  // Keep adding
1471  averaged_recovered_vort[nod_pt]+=recovered_vort;
1472  count[nod_pt]++;
1473  }
1474  }
1475 
1476  // Cleanup
1477  delete recovered_vorticity_coefficient_pt;
1478  recovered_vorticity_coefficient_pt=0;
1479  }
1480 
1481  //Loop over all nodes to actually work out the average
1482  unsigned nnod=mesh_pt->nnode();
1483  for(unsigned j=0;j<nnod;j++)
1484  {
1485  Node* nod_pt=mesh_pt->node_pt(j);
1486 
1487  //Calculate the values of the smoothed vorticity
1488  averaged_recovered_vort[nod_pt]/=count[nod_pt];
1489 
1490  //Assign smoothed vorticity to nodal values
1491  nod_pt->set_value(Smoothed_vorticity_index+deriv,
1492  averaged_recovered_vort[nod_pt]);
1493  }
1494 
1495 
1496  // Start again
1497  count.clear();
1498 
1499 
1500  } // end of loop over derivatives
1501 
1502  // Cleanup
1503  for (typename std::map<Node*,Vector<ELEMENT*>*>::iterator it=
1504  adjacent_elements_pt.begin();
1505  it!=adjacent_elements_pt.end();it++)
1506  {
1507  delete it->second;
1508  }
1509 
1510  oomph_info << "Time for vorticity recovery: "
1511  << TimingHelpers::timer()-t_start
1512  << " sec " << std::endl;
1513  }
1514 
1515  private:
1516 
1517  /// Order of recovery polynomials
1518  unsigned Recovery_order;
1519 
1520 };
1521 
1522 #endif
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
double vorticity_error_squared(const unsigned &i)
Compute the element's contribution to the (squared) L2 norm of the difference between exact and smoot...
void shape_rec(const Vector< double > &x, Vector< double > &psi_r)
cstr elem_len * i
Definition: cfortran.h:607
void pin_smoothed_vorticity()
Pin all smoothed vorticity quantities.
unsigned nrecovery_order() const
ExactVorticityFctPt Exact_vorticity_fct_pt
void get_interpolated_values(const unsigned &t, const Vector< double > &s, Vector< double > &values)
Get the function value u in Vector. Note: Given the generality of the interface (this function is usu...
char t
Definition: cfortran.h:572
void get_raw_velocity_deriv(const Vector< double > &s, Vector< double > &dveloc_dx) const
Get raw derivative of velocity.
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
void get_raw_vorticity_third_deriv(const Vector< double > &s, Vector< double > &dvorticity_dxdy) const
OomphInfo oomph_info
void setup_patches(Mesh *&mesh_pt, std::map< Node *, Vector< ELEMENT * > * > &adjacent_elements_pt, Vector< Node * > &vertex_node_pt)
Setup patches: For each vertex node pointed to by nod_pt, adjacent_elements_pt[nod_pt] contains the p...
e
Definition: cfortran.h:575
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:383
unsigned & recovery_order()
Access function for order of recovery polynomials.
VorticitySmoother(const VorticitySmoother &)
Broken copy constructor.
ExactVorticityFctPt & exact_vorticity_fct_pt()
unsigned smoothed_vorticity_index() const
Index of smoothed vorticity – followed by derivatives.
virtual ~VorticitySmoother()
Empty virtual destructor.
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:995
void output(std::ostream &outfile, const unsigned &nplot)
Overloaded output fct: Output veloc, pressure, smoothed vorticity.
void recover_vorticity(Mesh *mesh_pt)
Recover vorticity from patches.
void(* ExactVorticityFctPt)(const Vector< double > &x, double &vort, Vector< double > &dvort_dx, Vector< double > &dvort_dxdy, Vector< double > &dvort_dxdxdy, Vector< double > &dveloc_dx)
void vorticity_and_its_derivs(const Vector< double > &s, Vector< double > &veloc, double &vort, Vector< double > &dvort_dx, Vector< double > &dvort_dxdy, Vector< double > &dvort_dxdxdy, Vector< double > &dveloc_dx)
Compute smoothed vorticity and its derivatives.
void recover_vorticity(Mesh *mesh_pt, DocInfo &doc_info)
static char t char * s
Definition: cfortran.h:572
unsigned Recovery_order
Order of recovery polynomials.
unsigned required_nvalue(const unsigned &n) const
Number of values required at local node n. In order to simplify matters, we allocate storage for pres...
void output_smoothed_vorticity(std::ostream &outfile, const unsigned &nplot)
Output veloc, smoothed vorticity and derivatives.
ExactVorticityFctPt exact_vorticity_fct_pt() const
double timer()
returns the time in seconds after some point in past
void shape(const double &s, double *Psi)
Definition for 1D Lagrange shape functions. The value of all the shape functions at the local coordin...
Definition: shape.h:549
void get_interpolated_values(const Vector< double > &s, Vector< double > &values)
Get the function value u in Vector. Note: Given the generality of the interface (this function is usu...
unsigned ncont_interpolated_values() const
Number of continuously interpolated values:
void broken_assign(const std::string &class_name)
Issue error message and terminate execution.
VorticitySmoother(const unsigned &recovery_order)
Constructor: Set order of recovery shape functions.
Integral * integral_rec(const bool &is_q_mesh)
unsigned Number_of_values_per_field
Number of values per field: u,v,p,omega,d/dx,d/dy, d^2/dx^2,d^2/dxdy, d^2/dy^2,.
void get_recovered_vorticity_in_patch(const Vector< ELEMENT * > &patch_el_pt, const unsigned &num_recovery_terms, Vector< double > *&recovered_vorticity_coefficient_pt, unsigned &n_deriv)
void output_analytical_veloc_and_vorticity(std::ostream &outfile, const unsigned &nplot)
Output exact veloc, vorticity, derivs and indicator based on functions specified by two function poin...
void operator=(const VorticitySmoother &)
Broken assignment operator.
void get_raw_vorticity_second_deriv(const Vector< double > &s, Vector< double > &dvorticity_dxdy) const
Get raw derivative of smoothed derivative vorticity.
void get_raw_vorticity_deriv(const Vector< double > &s, Vector< double > &dvorticity_dx) const
Get raw derivative of smoothed vorticity.
Smoother for vorticity in 2D.
unsigned Smoothed_vorticity_index
Index of smoothed vorticity – followed by derivatives.