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// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26// Header file for Navier Stokes elements
27
28#ifndef OOMPH_VORTICITY_SMOOTHER_HEADER
29#define OOMPH_VORTICITY_SMOOTHER_HEADER
30
31// Config header generated by autoconfig
32#ifdef HAVE_CONFIG_H
33#include <oomph-lib-config.h>
34#endif
35
36
37//===============================================================
38// WARNING: THIS IS WORK IN PROGRESS -- ONLY USED IN 2D SO FAR
39//===============================================================
40namespace oomph
41{
42 //===========================================================================
43 /// Namespace with helper functions for (2D) vorticity (and derivatives)
44 /// recovery
45 //===========================================================================
46 namespace VorticityRecoveryHelpers
47 {
48 //========================================================================
49 /// Class to indicate which derivatives of the vorticity/
50 /// velocity we want to recover. We choose to immediately instantiate
51 /// an object of this class by dropping the semi-colon after the class
52 /// description.
53 //========================================================================
55 {
56 public:
57 /// Constructor
59 {
60 // Set the default (max) order of vorticity derivative to recover
62
63 // Set the default (max) order of velocity derivative to recover
65 }
66
67 /// The maximum order of derivatives calculated in the vorticity recovery
69 {
70 // Return the appropriate value
72 } // End of get_maximum_order_of_vorticity_derivative
73
74 /// The maximum order of derivatives calculated in the velocity recovery
76 {
77 // Return the appropriate value
79 } // End of get_maximum_order_of_velocity_derivative
80
81 /// The maximum order of derivatives calculated in the vorticity recovery
83 {
84 // Make sure the user has supplied a valid input value
85 if ((max_deriv < -1) || (max_deriv > 3))
86 {
87 // Throw an error
88 throw OomphLibError(
89 "Invalid input value! Should be between -1 and 3!",
90 OOMPH_CURRENT_FUNCTION,
91 OOMPH_EXCEPTION_LOCATION);
92 }
93
94 // Return the appropriate value
96
97 // Calculate the value of Number_of_values_per_field with the updated
98 // value of Maximum_order_of_vorticity_derivative
100 } // End of set_maximum_order_of_vorticity_derivative
101
102 /// The maximum order of derivatives calculated in the velocity recovery
104 {
105 // Make sure the user has supplied a valid input value. Note, unlike the
106 // vorticity, we always output the zeroth derivative of the velocity
107 // so we don't use -1 as an input
108 if ((max_deriv < 0) || (max_deriv > 1))
109 {
110 // Throw an error
111 throw OomphLibError("Invalid input value! Should be between 0 and 3!",
112 OOMPH_CURRENT_FUNCTION,
113 OOMPH_EXCEPTION_LOCATION);
114 }
115
116 // Return the appropriate value
118
119 // Calculate the value of Number_of_values_per_field with the updated
120 // value of Maximum_order_of_velocity_derivative
122 } // End of set_maximum_order_of_vorticity_derivative
123
124 /// Calculates the number of values per field given the number of
125 /// vorticity and velocity derivatives to recover (stored as private data)
127 {
128 // Output: u,v,p
130
131 // Loop over the vorticity derivatives
132 for (unsigned i = 0;
134 i++)
135 {
136 // Update the number of values per field
138 }
139
140 // Loop over the velocity derivatives
141 for (unsigned i = 1;
143 i++)
144 {
145 // Update the number of values per field
147 }
148 } // End of calculate_number_of_values_per_field
149
150 /// Helper function that determines the number of n-th order
151 /// partial derivatives in d-dimensions. Specifically there are
152 /// (n+d-1)(choose)(d-1) possible n-th order partial derivatives in
153 /// d-dimensions. Implementation makes use of the code found at:
154 /// www.geeksforgeeks.org/space-and-time-efficient-binomial-coefficient/
155 unsigned npartial_derivative(const unsigned& n) const
156 {
157 // This will only work in 2D so n_dim is always 2
158 unsigned n_dim = 2;
159
160 // Calculate m
161 unsigned n_bins = n + n_dim - 1;
162
163 // Calculate k
164 unsigned k = n_dim - 1;
165
166 // Initialise the result
167 unsigned value = 1;
168
169 // Since C(n_bins,k)=C(n_bins,n_bins-k)
170 if (k > n_bins - k)
171 {
172 // Replace k
173 k = n_bins - k;
174 }
175
176 // Calculate [n_bins*(n_bins-1)*...*(n_bins-k+1)]/[k*(k-1)*...*1]
177 for (unsigned i = 0; i < k; ++i)
178 {
179 // First update
180 value *= (n_bins - i);
181
182 // Second update
183 value /= (i + 1);
184 }
185
186 // Return the result
187 return value;
188 } // End of npartial_derivative
189
190 /// Number of continuously interpolated values:
192 {
193 // Return the number of values used per field
195 } // End of ncont_interpolated_values
196
197 private:
198 /// Number of values per field; how many of the following do we
199 /// want: u,v,p,omega,d/dx,d/dy, d^2/dx^2,d^2/dxdy,d^2/dy^2,
200 /// d^3/dx^3,d^3/dx^2dy,d^3/dxdy^2,d^3/dy^3,
201 /// du/dx,du/dy,dv/dx,dv/dy
203
204 /// Maximum number of derivatives to retain in the vorticity
205 /// recovery. Note, the value -1 means we ONLY output u,v[,w],p.
207
208 /// Maximum number of derivatives to retain in the velocity
209 /// recovery. Note, the value 0 means we don't calculate the
210 /// derivatives of the velocity
213
214
215 } // namespace VorticityRecoveryHelpers
216
217
218 //===============================================
219 /// Overloaded element that allows projection of
220 /// vorticity.
221 //===============================================
222 template<class ELEMENT>
223 class VorticitySmootherElement : public virtual ELEMENT
224 {
225 public:
226 /// Constructor
228 {
229 // Only done for 2D (3D is far more difficult -- recovering a vector field
230 // instead of a scalar field)
231 N_dim = 2;
232
233 // Index of smoothed vorticity ([0]:u, [1]:v, [2]:p ==> [3]:w)
235
236 // The current maximum order of vorticity derivatives we can recover. This
237 // is currently 10 as we can recover from the zeroth to third derivative
239
240 // The current maximum order of velocity derivatives we can recover.
241 // Currently this is 4 as we only recover first derivatives.
243
244 // Recover as many
248
249 // The default is to recover the velocity derivatives
253
254 // Set the number of values per field
257
258 // Pointer to fct that specifies exact vorticity and
259 // derivs (for validation).
261 }
262
263 /// Typedef for pointer to function that specifies the exact vorticity
264 /// and derivatives (for validation)
265 typedef void (*ExactVorticityFctPt)(
266 const Vector<double>& x, Vector<Vector<double>>& vort_and_derivs);
267
268 /// Helper function to create a container for the vorticity and
269 /// its partial derivatives. If the user wishes to output everything then
270 /// this also creates space for the velocity derivatives too. This function
271 /// has been written to allow for generalisation of the storage without
272 /// (hopefully!) affecting too much
274 const
275 {
276 // Maximum number of vorticity derivatives
277 unsigned n_vort_deriv = Maximum_order_of_vorticity_derivative;
278
279 // Maximum number of velocity derivatives
280 unsigned n_veloc_deriv = Maximum_order_of_velocity_derivative;
281
282 // The number of vectors in the output
283 unsigned n_vector = n_vort_deriv + n_veloc_deriv + 1;
284
285 // Container for the vorticity and its derivatives
286 Vector<Vector<double>> vort_and_derivs(n_vector);
287
288 // Loop over the entries of the vector
289 for (unsigned i = 0; i < n_vort_deriv + 1; i++)
290 {
291 // Get the number of partial derivatives of the vorticity
292 vort_and_derivs[i].resize(npartial_derivative(i), 0.0);
293 }
294
295 // Loop over the entries of the vector
296 for (unsigned i = n_vort_deriv + 1; i < n_vort_deriv + n_veloc_deriv + 1;
297 i++)
298 {
299 // Resize the container and initialise all entries to zero
300 vort_and_derivs[i].resize(2 * npartial_derivative(i - n_vort_deriv),
301 0.0);
302 }
303
304 // Return the container
305 return vort_and_derivs;
306 } // End of create_container_for_vorticity_and_derivatives
307
308 /// Helper function that, given the local dof number of the i-th
309 /// vorticity or velocity derivative, returns the index in the container
310 /// that stores the corresponding value.
311 /// Note 1: this function goes hand-in-hand with
312 /// create_container_for_vorticity_and_derivatives()
313 /// Note 2: i=0: vorticity itself;
314 /// i>0: vorticity derivatives
315 std::pair<unsigned, unsigned> vorticity_dof_to_container_id(
316 const unsigned& i) const
317 {
318 // Maximum number of vorticity derivatives (that we actually want)
319 unsigned n_vort_deriv = Maximum_order_of_vorticity_derivative;
320
321 // Maximum number of velocity derivatives (that we actually want)
322 unsigned n_veloc_deriv = Maximum_order_of_velocity_derivative;
323
324#ifdef PARANOID
325 // We cannot calculate this if we're not recovering anything
326 if (n_vort_deriv + n_veloc_deriv + 1 == 0)
327 {
328 // Throw an error
329 throw OomphLibError(
330 "Not recovering anything so this shouldn't be called.",
331 OOMPH_CURRENT_FUNCTION,
332 OOMPH_EXCEPTION_LOCATION);
333 }
334#endif
335
336 // Maximum order of vorticity derivative (that can be recovered)
337 unsigned max_vort_order =
339
340 // Maximum order of velocity derivative (that can be recovered)
341 unsigned max_veloc_order =
343
344 // Maximum number of recoverable vorticity terms
345 unsigned max_vort_recov = 0;
346
347 // Maximum number of recoverable velocity terms
348 unsigned max_veloc_recov = 0;
349
350 // Loop over the entries of the vector
351 for (unsigned j = 0; j < max_vort_order + 1; j++)
352 {
353 // Get the number of partial derivatives of the vorticity
354 max_vort_recov += npartial_derivative(j);
355 }
356
357 // Loop over the entries of the vector
358 for (unsigned j = 1; j < max_veloc_order + 1; j++)
359 {
360 // Get the number of partial derivatives of the vorticity
361 max_veloc_recov += 2 * npartial_derivative(j);
362 }
363
364 // Create a pair to store the output
365 std::pair<unsigned, unsigned> container_id;
366
367 // Store the number of derivatives we've gone over
368 unsigned nprev_deriv = 0;
369
370 // The number of derivatives available of a given order
371 unsigned n_deriv = 0;
372
373 // If the dof number i is associated with a vorticity derivative
374 if (i < max_vort_recov)
375 {
376 // Loop over the vorticity vectors
377 for (unsigned jj = 0; jj < n_vort_deriv + 1; jj++)
378 {
379 // Number of derivatives of order jj
380 n_deriv = npartial_derivative(jj);
381
382 // Update nprev_derivs
383 nprev_deriv += n_deriv;
384
385 // If this is true then we need the jj-th vector
386 if (i < nprev_deriv)
387 {
388 // Assign the first index
389 container_id.first = jj;
390
391 // Index of the entry
392 container_id.second = i - (nprev_deriv - n_deriv);
393
394 // We're done
395 return container_id;
396 }
397 } // for (unsigned jj=0;jj<n_vort_deriv+1;jj++)
398 }
399 // If the dof number i is associated with a velocity derivative
400 else if (i < max_vort_recov + max_veloc_recov)
401 {
402 // Initialise the value of nprev_deriv (depending on how many
403 // vorticity derivatives we can recover)
404 nprev_deriv = max_vort_recov;
405
406 // Loop over the velocity vectors
407 for (unsigned jj = n_vort_deriv + 1;
408 jj < n_vort_deriv + n_veloc_deriv + 1;
409 jj++)
410 {
411 // Number of velocity derivatives
412 n_deriv = 2 * npartial_derivative(jj - n_vort_deriv);
413
414 // Update nprev_derivs
415 nprev_deriv += n_deriv;
416
417 // If this is true then we need the jj-th vector
418 if (i < nprev_deriv)
419 {
420 // Assign the first index
421 container_id.first = jj;
422
423 // Index of the entry
424 container_id.second = i - (nprev_deriv - n_deriv);
425
426 // We're done
427 return container_id;
428 }
429 } // for (unsigned jj=n_vort_deriv+1;jj<n_veloc_deriv;jj++)
430 }
431 else
432 {
433 // Create an ostringstream object to create an error message
434 std::ostringstream error_message_stream;
435
436 // Create an error message
437 error_message_stream << "Dof number " << i << " not associated "
438 << "with a vorticity or velocity derivative!"
439 << std::endl;
440
441 // Throw an error
442 throw OomphLibError(error_message_stream.str(),
443 OOMPH_CURRENT_FUNCTION,
444 OOMPH_EXCEPTION_LOCATION);
445 } // if (i<max_vort_recov)
446
447 // We'll never get here but need to return something
448 return container_id;
449 } // End of vorticity_dof_to_container_id
450
451
452 /// Helper function that, given the local dof number of the i-th
453 /// vorticity or velocity derivative, returns the index in the container
454 /// that stores the corresponding value.
455 /// Note 1: this function goes hand-in-hand with
456 /// create_container_for_vorticity_and_derivatives()
457 /// Note 2: The input to this function is the i associated with the STORED
458 /// nodal dof value. For example, if we're only recovering the
459 /// velocity derivatives then i=0 is associated with du/dx
460 std::pair<unsigned, unsigned> recovered_dof_to_container_id(
461 const unsigned& i) const
462 {
463 // Create a pair to store the output
464 std::pair<unsigned, unsigned> container_id;
465
466 // Find which derivative to calculate
467 unsigned derivative_index = i;
468
469 // Variable to store the index of the vector (initialise the value to -1
470 // so we know whether or not the value has been set)
471 int vector_index = -1;
472
473 // Get the number of vorticity derivatives to recover
474 unsigned n_vort_derivs = Maximum_order_of_vorticity_derivative;
475
476 // Get the number of vorticity derivatives to recover
477 unsigned n_veloc_derivs = Maximum_order_of_velocity_derivative;
478
479 // Loop over the entries of the vector
480 for (unsigned i = 0; i < n_vort_derivs + 1; i++)
481 {
482 // Get the number of partial derivatives of the vorticity
483 if (derivative_index < npartial_derivative(i))
484 {
485 // We're on the i-th vector
486 vector_index = i;
487
488 // We're done
489 break;
490 }
491 else
492 {
493 // Decrease the value of derivative_index
494 derivative_index -= npartial_derivative(i);
495 }
496 } // for (unsigned i=0;i<n_vort_derivs+1;i++)
497
498 // If the vector_index variable value hasn't been found yet
499 if (vector_index == -1)
500 {
501 // Loop over the entries of the vector
502 for (unsigned i = 1; i < n_veloc_derivs + 1; i++)
503 {
504 // Get the number of partial derivatives of the vorticity
505 if (derivative_index < 2 * npartial_derivative(i))
506 {
507 // We're on the i-th vector
508 vector_index = n_vort_derivs + i;
509
510 // We're done
511 break;
512 }
513 else
514 {
515 // Decrease the value of derivative_index
516 derivative_index -= 2 * npartial_derivative(i);
517 }
518 } // for (unsigned i=1;i<n_veloc_derivs+1;i++)
519 } // if (vector_index==-1)
520
521#ifdef PARANOID
522 // Sanity check: if the value hasn't been set yet, something's wrong
523 if (vector_index == -1)
524 {
525 // Create an ostringstream object to create an error message
526 std::ostringstream error_message_stream;
527
528 // Create an error message
529 error_message_stream << "Value of vector_index has not been set. "
530 << "Something's wrong!";
531
532 // Throw an error
533 throw OomphLibError(error_message_stream.str(),
534 OOMPH_CURRENT_FUNCTION,
535 OOMPH_EXCEPTION_LOCATION);
536 }
537#endif
538
539 // Assign the first entry of container_id
540 container_id.first = vector_index;
541
542 // Assign the second entry of container_id
543 container_id.second = derivative_index;
544
545 // We'll never get here but need to return something
546 return container_id;
547 } // End of recovered_dof_to_container_id
548
549 /// Given the STORED dof number, this function returns the global
550 /// recovered number. For example, if we only want to recover the velocity
551 /// derivatives then the stored dof number of du/dy is 1 (as 0 is associated
552 /// with du/dx). The global recovered number is 11 (as there are currently
553 /// 10 vorticity derivatives that can be recovered and du/dy is the second
554 /// velocity derivative we can recover).
555 unsigned stored_dof_to_recoverable_dof(const unsigned& i) const
556 {
557 // Get the ID in the storage associated with i-th recovered dof
558 std::pair<unsigned, unsigned> id =
560
561 // Vector entry index
562 unsigned vector_index = id.first;
563
564 // Derivative index
565 unsigned deriv_index = id.second;
566
567 // The index we want
568 unsigned index = 0;
569
570 // Maximum possible order of vorticity derivatives that can be recovered
571 unsigned max_vort_recov =
573
574 // If we're dealing with a vorticity derivative
575 if (vector_index < max_vort_recov + 1)
576 {
577 // Holds the number of derivatives of lower order
578 unsigned n_prev_deriv = 0;
579
580 // Loop over the derivative orders
581 for (unsigned j = 0; j < vector_index; j++)
582 {
583 // Increment n_prev_deriv
584 n_prev_deriv += npartial_derivative(j);
585 }
586
587 // Update the value of index
588 index += n_prev_deriv + deriv_index;
589 }
590 // We're dealing with derivatives of the velocity
591 else
592 {
593 // Holds the number of vorticity derivatives
594 unsigned n_prev_deriv = 0;
595
596 // Loop over the derivative orders
597 for (unsigned j = 0; j < max_vort_recov + 1; j++)
598 {
599 // Increment n_prev_deriv
600 n_prev_deriv += npartial_derivative(j);
601 }
602
603 // What is the user-chosen maximum vorticity derivative to recover?
604 unsigned n_vort_deriv = Maximum_order_of_vorticity_derivative;
605
606 // Loop over the derivative orders
607 for (unsigned j = 1; j < vector_index - n_vort_deriv; j++)
608 {
609 // Increment n_prev_deriv
610 n_prev_deriv += 2 * npartial_derivative(j);
611 }
612
613 // Update the value of index
614 index += n_prev_deriv + deriv_index;
615 } // if (vector_index<max_vort_recov)
616
617 // Return the value of index
618 return index;
619 } // End of stored_dof_to_recoverable_dof
620
621 /// The maximum order of vorticity derivative that can be recovered.
622 /// This is set in the constructor and should NOT be changed during the
623 /// running of the code. As such, a set...() function is not provided.
624 /// DRAIG: Leave get_ prefix?
626 {
627 // Return the appropriate value
629 } // End of get_maximum_order_of_recoverable_vorticity_derivative
630
631 /// The maximum order of velocity derivative that can be recovered.
632 /// This is set in the constructor and should NOT be changed during the
633 /// running of the code. As such, a set...() function is not provided.
634 /// DRAIG: Leave get_ prefix?
636 {
637 // Return the appropriate value
639 } // End of get_maximum_order_of_recoverable_velocity_derivative
640
641 /// The maximum order of derivatives calculated in the vorticity
642 /// recovery. Note, this value can only be set through the namespace
643 /// VorticityRecoveryHelpers.
644 /// DRAIG: Leave get_ prefix?
646 {
647 // Return the appropriate value
649 } // End of get_maximum_order_of_vorticity_derivative
650
651 /// The maximum order of derivatives calculated in the velocity
652 /// recovery. Note, this value can only be set through the namespace
653 /// VorticityRecoveryHelpers.
654 /// DRAIG: Leave get_ prefix?
656 {
657 // Return the appropriate value
659 } // End of get_maximum_order_of_velocity_derivative
660
661 /// The number of terms calculated in the vorticity recovery. Also
662 /// includes the zeroth derivative, i.e. the vorticity itself
664 {
665 // The number of terms recovered
666 unsigned n_terms = 0;
667
668 // Loop over the derivatives
669 for (unsigned i = 0;
671 i++)
672 {
673 // Update n_terms
674 n_terms += npartial_derivative(i);
675 }
676
677 // Return the appropriate value
678 return n_terms;
679 } // End of nvorticity_derivatives_to_recover
680
681 /// The number of derivatives calculated in the velocity recovery. This does
682 /// NOT include the zeroth derivatives as they are not "recovered"
684 {
685 // The number of terms recovered
686 unsigned n_terms = 0;
687
688 // Loop over the derivatives
689 for (unsigned i = 1;
691 i++)
692 {
693 // Update n_terms
694 n_terms += 2 * npartial_derivative(i);
695 }
696
697 // Return the appropriate value
698 return n_terms;
699 } // End of nvelocity_derivatives_to_recover
700
701 /// Call the function written in VorticityRecoveryHelpers
702 unsigned npartial_derivative(const unsigned& n) const
703 {
704 // Return the result
706 } // End of npartial_derivative
707
708 /// Access function: Pointer to function that specifies exact
709 /// vorticity and derivatives (for validation).
711 {
712 // Return the address of the function pointer
714 } // End of exact_vorticity_fct_pt
715
716 /// Access function: Pointer to function that specifies exact
717 /// vorticity and derivatives (for validation) -- const version
719 {
720 // Return the address of the function pointer
722 } // End of exact_vorticity_fct_pt
723
724 /// Index of smoothed vorticity -- followed by derivatives
726 {
727 // Return the value of the Smoothed_vorticity_index
729 } // End of smoothed_vorticity_index
730
731 /// Number of values required at local node n. In order to simplify
732 /// matters, we allocate storage for pressure variables at all the nodes
733 /// and then pin those that are not used.
734 unsigned required_nvalue(const unsigned& n) const
735 {
736 // Return the number of values used per field
738 } // End of required_nvalue
739
740 /// Number of continuously interpolated values:
742 {
743 // Return the number of values used per field
745 } // End of ncont_interpolated_values
746
747 /// Get the function value u in Vector.
748 /// Note: Given the generality of the interface (this function is usually
749 /// called from black-box documentation or interpolation routines), the
750 /// values Vector sets its own size in here.
752 Vector<double>& values)
753 {
754 // Get the value at the current time
755 unsigned t = 0;
756
757 // Get the interpolated values
758 get_interpolated_values(t, s, values);
759 } // End of get_interpolated_values
760
761 /// Get the function value u in Vector.
762 /// Note: Given the generality of the interface (this function is usually
763 /// called from black-box documentation or interpolation routines), the
764 /// values Vector sets its own size in here.
765 void get_interpolated_values(const unsigned& t,
766 const Vector<double>& s,
767 Vector<double>& values)
768 {
769 // Set size of vector and initialise all entries to zero
770 values.resize(Number_of_values_per_field, 0.0);
771
772 // Find out how many nodes there are
773 unsigned n_node = this->nnode();
774
775 // Shape functions
776 Shape psif(n_node);
777
778 // Get the value of the shape function at the local coordinate s
779 this->shape(s, psif);
780
781 // Calculate velocities: values[0],...
782 for (unsigned i = 0; i < N_dim; i++)
783 {
784 // Get the index at which the i-th velocity is stored
785 unsigned u_nodal_index = this->u_index_nst(i);
786
787 // Loop over the nodes
788 for (unsigned l = 0; l < n_node; l++)
789 {
790 // Update the i-th entry of the value vector
791 values[i] += this->nodal_value(t, l, u_nodal_index) * psif[l];
792 }
793 } // for (unsigned i=0;i<N_dim;i++)
794
795 // Calculate pressure: values[N_dim] (no history is carried in the
796 // pressure)
797 values[N_dim] = this->interpolated_p_nst(s);
798 } // End of get_interpolated_values
799
800 /// Pin all smoothed vorticity quantities
802 {
803 // Get the number of nodes in the element
804 unsigned nnod = this->nnode();
805
806 // Loop over the nodes
807 for (unsigned j = 0; j < nnod; j++)
808 {
809 // Make a pointer to the j-th node in the element
810 Node* nod_pt = this->node_pt(j);
811
812 // Loop over the fields
813 for (unsigned i = Smoothed_vorticity_index;
815 i++)
816 {
817 // Pin the i-th field at this node
818 nod_pt->pin(i);
819 }
820 } // for (unsigned j=0;j<nnod;j++)
821 } // End of pin_smoothed_vorticity
822
823 /// Output exact velocity, vorticity, derivatives and indicator
824 /// based on functions specified by two function pointers
825 void output_analytical_veloc_and_vorticity(std::ostream& outfile,
826 const unsigned& nplot)
827 {
828 // Vector of local coordinates
829 Vector<double> s(N_dim, 0.0);
830
831 // Global coordinates container
833
834 // Get the number of nodes in this element
835 unsigned n_node = this->nnode();
836
837 // Shape functions
838 Shape psif(n_node);
839
840 // Tecplot header info
841 outfile << this->tecplot_zone_string(nplot);
842
843 // Get the number of plot points
844 unsigned num_plot_points = this->nplot_points(nplot);
845
846 // Create a container for the vorticity and its derivatives
847 Vector<Vector<double>> vort_and_derivs =
849
850 // Loop over plot points
851 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
852 {
853 // Get local coordinates of plot point
854 this->get_s_plot(iplot, nplot, s);
855
856 // Loop over the coordinates
857 for (unsigned i = 0; i < N_dim; i++)
858 {
859 // Get the interpolated coordinate value at this local coordinate
860 x[i] = this->interpolated_x(s, i);
861
862 // Output the i-th coordinate value to file
863 outfile << x[i] << " ";
864 }
865
866 // Fake velocity and pressure
867 outfile << "0.0 0.0 0.0 ";
868
869 // Get the vorticity and its derivatives (and the derivatives of the
870 // velocity if required)
871 Exact_vorticity_fct_pt(x, vort_and_derivs);
872
873 // Number of vectors
874 unsigned n_vector = vort_and_derivs.size();
875
876 // Loop over the number of vectors we're outputting from
877 for (unsigned i = 0; i < n_vector; i++)
878 {
879 // The number of entries in the vector
880 unsigned i_entries = vort_and_derivs[i].size();
881
882 // Loop over the entries in the vector
883 for (unsigned j = 0; j < i_entries; j++)
884 {
885 // Output the smoothed vorticity to file
886 outfile << (vort_and_derivs[i])[j] << " ";
887 }
888 } // for (unsigned i=0;i<n_vector;i++)
889
890 // Finish the line
891 outfile << std::endl;
892 } // for (unsigned iplot=0;iplot<num_plot_points;iplot++)
893
894 // Write tecplot footer (e.g. FE connectivity lists)
895 this->write_tecplot_zone_footer(outfile, nplot);
896 } // End of output_analytical_veloc_and_vorticity
897
898 /// Output the velocity, smoothed vorticity and derivatives
899 void output_smoothed_vorticity(std::ostream& outfile, const unsigned& nplot)
900 {
901 // Vector of local coordinates
902 Vector<double> s(N_dim, 0.0);
903
904 // Vector of global coordinates
905 Vector<double> x(N_dim, 0.0);
906
907 // Get the number of nodes in the element
908 unsigned n_node = this->nnode();
909
910 // Shape functions
911 Shape psif(n_node);
912
913 // Tecplot header info
914 outfile << this->tecplot_zone_string(nplot);
915
916 // Create a container for the vorticity and its derivatives
917 Vector<Vector<double>> vort_and_derivs =
919
920 // Vector of velocities
921 Vector<double> velocity(N_dim, 0.0);
922
923 // Get the number of plot points
924 unsigned num_plot_points = this->nplot_points(nplot);
925
926 // Loop over plot points
927 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
928 {
929 // Get local coordinates of plot point
930 this->get_s_plot(iplot, nplot, s);
931
932 // Get vorticity and its derivatives (reconstructed)
933 vorticity_and_its_derivs(s, vort_and_derivs);
934
935 // Loop over the coordinates
936 for (unsigned i = 0; i < N_dim; i++)
937 {
938 // Get the i-th interpolated coordinate at a given local coordinate
939 x[i] = this->interpolated_x(s, i);
940
941 // Output the interpolated coordinate to file
942 outfile << x[i] << " ";
943 }
944
945 // Loop over the coordinates
946 for (unsigned i = 0; i < N_dim; i++)
947 {
948 // Output the i-th velocity component
949 outfile << this->interpolated_u_nst(s, i) << " ";
950 }
951
952 // Pressure
953 outfile << this->interpolated_p_nst(s) << " ";
954
955 // Output the smoothed vorticity derivatives
956 // (d/dx, d/dy, d^2/dx^2, d^2/dxdy, d^2/dy^2
957 // d^3/dx^3, d^3/dx^2dy, d^3/dxdy^2, d^3/dy^3
958 // du/dx, du/dy, dv/dx, dv/dy):
959
960 // Number of vectors
961 unsigned n_vector = vort_and_derivs.size();
962
963 // Loop over the number of vectors we're outputting from
964 for (unsigned i = 0; i < n_vector; i++)
965 {
966 // The number of entries in the vector
967 unsigned i_entries = vort_and_derivs[i].size();
968
969 // Loop over the entries in the vector
970 for (unsigned j = 0; j < i_entries; j++)
971 {
972 // Output the smoothed vorticity to file
973 outfile << (vort_and_derivs[i])[j] << " ";
974 }
975 } // for (unsigned i=0;i<n_vector;i++)
976
977 // End the line
978 outfile << std::endl;
979 } // for (unsigned iplot=0;iplot<num_plot_points;iplot++)
980
981 // Write tecplot footer (e.g. FE connectivity lists)
982 this->write_tecplot_zone_footer(outfile, nplot);
983 } // End of output_smoothed_vorticity
984
985 /// Number of scalars/fields output by this element. Re-implements
986 /// broken virtual function in base class.
987 unsigned nscalar_paraview() const
988 {
989 // Return the number of continuously interpolated values
991 } // End of nscalar_paraview
992
993 /// Write values of the i-th scalar field at the plot points. Needs
994 /// to be implemented for each new specific element type.
995 void scalar_value_paraview(std::ofstream& file_out,
996 const unsigned& i,
997 const unsigned& nplot) const
998 {
999 // Vector of local coordinates
1000 Vector<double> s(N_dim, 0.0);
1001
1002 // Get the number of plot points
1003 unsigned num_plot_points = this->nplot_points_paraview(nplot);
1004
1005 // Create a container for the vorticity and its derivatives
1006 Vector<Vector<double>> vort_and_derivs =
1008
1009 // Loop over plot points
1010 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
1011 {
1012 // Get local coordinates of plot point
1013 this->get_s_plot(iplot, nplot, s);
1014
1015 // Velocities
1016 if (i < N_dim)
1017 {
1018 // Output the interpolated velocity component
1019 file_out << this->interpolated_u_nst(s, i) << std::endl;
1020 }
1021 // Pressure
1022 else if (i == N_dim)
1023 {
1024 // Output the interpolated pressure value
1025 file_out << this->interpolated_p_nst(s) << std::endl;
1026 }
1027 // Output the vorticity and required derivatives
1028 else if (i < nscalar_paraview())
1029 {
1030 // Get vorticity and its derivatives (reconstructed)
1031 vorticity_and_its_derivs(s, vort_and_derivs);
1032
1033 // Get the ID in the storage associated with i-th recovered dof
1034 std::pair<unsigned, unsigned> id =
1036
1037 // Output the appropriate value
1038 file_out << (vort_and_derivs[id.first])[id.second] << std::endl;
1039 }
1040 // Never get here
1041 else
1042 {
1043#ifdef PARANOID
1044 // Using a std::stringstream object to create a string
1045 std::stringstream error_stream;
1046
1047 // Create the error message
1048 error_stream << "These VorticitySmoother elements only store "
1049 << ncont_interpolated_values() << " fields, "
1050 << "but i is currently: " << i << std::endl;
1051
1052 // Throw an error
1053 throw OomphLibError(error_stream.str(),
1054 OOMPH_CURRENT_FUNCTION,
1055 OOMPH_EXCEPTION_LOCATION);
1056#endif
1057 }
1058 } // for (unsigned iplot=0;iplot<num_plot_points;iplot++)
1059 } // End of scalar_value_paraview
1060
1061 /// Name of the i-th scalar field. Default implementation
1062 /// returns V1 for the first one, V2 for the second etc. Can (should!) be
1063 /// overloaded with more meaningful names in specific elements.
1064 std::string scalar_name_paraview(const unsigned& i) const
1065 {
1066 // Velocities
1067 if (i < N_dim)
1068 {
1069 // Return the velocity name
1070 return "Velocity " + StringConversion::to_string(i);
1071 }
1072 // Pressure
1073 else if (i == N_dim)
1074 {
1075 // Return the appropriate string
1076 return "Pressure";
1077 }
1078 // Vorticity or its derivatives
1079 else if (i < nscalar_paraview())
1080 {
1081 // Calculate the appropriate value of index
1082 unsigned index =
1084
1085 // Switch statement is the easiest way to output smoothed vorticity
1086 // and velocity derivatives:
1087 // (d/dx, d/dy,
1088 // d^2/dx^2, d^2/dxdy, d^2/dy^2
1089 // d^3/dx^3, d^3/dx^2dy, d^3/dxdy^2, d^3/dy^3
1090 // du/dx, du/dy, dv/dx, dv/dy)
1091 switch (index)
1092 {
1093 case 3:
1094 return "w";
1095 break;
1096 case 4:
1097 return "dw/dx";
1098 break;
1099 case 5:
1100 return "dw/dy";
1101 break;
1102 case 6:
1103 return "d^2w/dx^2";
1104 break;
1105 case 7:
1106 return "d^2w/dxdy";
1107 break;
1108 case 8:
1109 return "d^2w/dy^2";
1110 break;
1111 case 9:
1112 return "d^3/dx^3";
1113 break;
1114 case 10:
1115 return "d^3/dx^2dy";
1116 break;
1117 case 11:
1118 return "d^3/dxdy^2";
1119 break;
1120 case 12:
1121 return "d^3/dy^3";
1122 break;
1123 case 13:
1124 return "du/dx";
1125 break;
1126 case 14:
1127 return "du/dy";
1128 break;
1129 case 15:
1130 return "dv/dx";
1131 break;
1132 case 16:
1133 return "dv/dy";
1134 break;
1135 default:
1136 oomph_info << "Never get here" << std::endl;
1137 abort();
1138 break;
1139 }
1140 }
1141 // Never get here
1142 else
1143 {
1144 std::stringstream error_stream;
1145 error_stream << "These Navier Stokes elements only store "
1146 << nscalar_paraview() << " fields,\n"
1147 << "but i is currently: " << i << std::endl;
1148
1149 // Throw the error
1150 throw OomphLibError(
1151 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1152 // Dummy return
1153 return " ";
1154 }
1155 } // End of scalar_name_paraview
1156
1157 /// Overloaded output function: Output velocity, pressure and the
1158 /// smoothed vorticity
1159 void output(std::ostream& outfile, const unsigned& nplot)
1160 {
1161 // Vector of local coordinates
1162 Vector<double> s(N_dim, 0.0);
1163
1164 // Global coordinates container
1165 Vector<double> x(N_dim, 0.0);
1166
1167 // Get the number of nodes in this element
1168 unsigned n_node = this->nnode();
1169
1170 // Shape functions
1171 Shape psif(n_node);
1172
1173 // Tecplot header info
1174 outfile << this->tecplot_zone_string(nplot);
1175
1176 // Get the number of plot points
1177 unsigned num_plot_points = this->nplot_points(nplot);
1178
1179 // Create a container for the vorticity and its derivatives
1180 Vector<Vector<double>> vort_and_derivs =
1182
1183 // Loop over plot points
1184 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
1185 {
1186 // Get local coordinates of plot point
1187 this->get_s_plot(iplot, nplot, s);
1188
1189 // Get shape fct
1190 this->shape(s, psif);
1191
1192 // Loop over the coordinates
1193 for (unsigned i = 0; i < N_dim; i++)
1194 {
1195 // Get the interpolated coordinate value at this local coordinate
1196 x[i] = this->interpolated_x(s, i);
1197
1198 // Output the i-th coordinate value to file
1199 outfile << x[i] << " ";
1200 }
1201
1202 // Loop over the coordinates
1203 for (unsigned i = 0; i < N_dim; i++)
1204 {
1205 // Output the i-th velocity component
1206 outfile << this->interpolated_u_nst(s, i) << " ";
1207 }
1208
1209 // Pressure
1210 outfile << this->interpolated_p_nst(s) << " ";
1211
1212 // Get vorticity and its derivatives (reconstructed)
1213 vorticity_and_its_derivs(s, vort_and_derivs);
1214
1215 // Number of vectors
1216 unsigned n_vector = vort_and_derivs.size();
1217
1218 // Loop over the number of vectors we're outputting from
1219 for (unsigned i = 0; i < n_vector; i++)
1220 {
1221 // The number of entries in the vector
1222 unsigned i_entries = vort_and_derivs[i].size();
1223
1224 // Loop over the entries in the vector
1225 for (unsigned j = 0; j < i_entries; j++)
1226 {
1227 // Output the smoothed vorticity to file
1228 outfile << (vort_and_derivs[i])[j] << " ";
1229 }
1230 } // for (unsigned i=0;i<n_vector;i++)
1231
1232 // Finish the line off
1233 outfile << std::endl;
1234 }
1235
1236 // Write tecplot footer (e.g. FE connectivity lists)
1237 this->write_tecplot_zone_footer(outfile, nplot);
1238 } // End of output
1239
1240
1241 /// Get raw derivative of velocity
1243 Vector<double>& dveloc_dx) const
1244 {
1245 // Find out how many nodes there are
1246 unsigned n_node = this->nnode();
1247
1248 // Set up memory for the shape functions
1249 Shape psif(n_node);
1250
1251 // Set up memory for the shape function derivatives
1252 DShape dpsifdx(n_node, 2);
1253
1254 // Call the derivatives of the shape and test functions
1255 this->dshape_eulerian(s, psif, dpsifdx);
1256
1257 // Initialise all entries to zero
1258 dveloc_dx.initialise(0.0);
1259
1260 // Loop over nodes
1261 for (unsigned l = 0; l < n_node; l++)
1262 {
1263 // Loop over derivative directions
1264 for (unsigned j = 0; j < 2; j++)
1265 {
1266 // Derivatives of the first velocity component
1267 dveloc_dx[j] += this->nodal_value(l, 0) * dpsifdx(l, j);
1268
1269 // Derivatives of the second velocity component
1270 dveloc_dx[j + 2] += this->nodal_value(l, 1) * dpsifdx(l, j);
1271 }
1272 } // for (unsigned l=0;l<n_node;l++)
1273 } // End of get_raw_velocity_deriv
1274
1275 /// Get raw derivative of velocity
1277 double& dveloc_dx,
1278 const unsigned& index) const
1279 {
1280 // Find out how many nodes there are
1281 unsigned n_node = this->nnode();
1282
1283 // Set up memory for the shape functions
1284 Shape psif(n_node);
1285
1286 // Set up memory for the shape function derivatives
1287 DShape dpsifdx(n_node, 2);
1288
1289 // Call the derivatives of the shape and test functions
1290 this->dshape_eulerian(s, psif, dpsifdx);
1291
1292 // Initialise value to zero
1293 dveloc_dx = 0.0;
1294
1295 // Use a switch statement
1296 switch (index)
1297 {
1298 case 0:
1299 // Loop over the nodes
1300 for (unsigned l = 0; l < n_node; l++)
1301 {
1302 // Derivatives of the first velocity component
1303 dveloc_dx += this->nodal_value(l, 0) * dpsifdx(l, 0);
1304 }
1305 break;
1306 case 1:
1307 // Loop over the nodes
1308 for (unsigned l = 0; l < n_node; l++)
1309 {
1310 // Derivatives of the first velocity component
1311 dveloc_dx += this->nodal_value(l, 0) * dpsifdx(l, 1);
1312 }
1313 break;
1314 case 2:
1315 // Loop over the nodes
1316 for (unsigned l = 0; l < n_node; l++)
1317 {
1318 // Derivatives of the second velocity component
1319 dveloc_dx += this->nodal_value(l, 1) * dpsifdx(l, 0);
1320 }
1321 break;
1322 case 3:
1323 // Loop over the nodes
1324 for (unsigned l = 0; l < n_node; l++)
1325 {
1326 // Derivatives of the second velocity component
1327 dveloc_dx += this->nodal_value(l, 1) * dpsifdx(l, 1);
1328 }
1329 break;
1330 default:
1331 oomph_info << "Never get here!" << std::endl;
1332 abort();
1333 }
1334 } // End of get_raw_velocity_deriv
1335
1336 /// Get raw derivative of smoothed vorticity
1338 Vector<double>& dvorticity_dx) const
1339 {
1340 // Find out how many nodes there are
1341 unsigned n_node = this->nnode();
1342
1343 // Set up memory for the shape functions
1344 Shape psif(n_node);
1345
1346 // Set up memory for the shape function derivatives
1347 DShape dpsifdx(n_node, 2);
1348
1349 // Call the derivatives of the shape and test functions
1350 this->dshape_eulerian(s, psif, dpsifdx);
1351
1352 // Initialise all entries to zero
1353 dvorticity_dx.initialise(0.0);
1354
1355 // Loop over nodes
1356 for (unsigned l = 0; l < n_node; l++)
1357 {
1358 // Loop over derivative directions
1359 for (unsigned j = 0; j < 2; j++)
1360 {
1361 // Calculate the x and y derivative
1362 dvorticity_dx[j] +=
1363 this->nodal_value(l, Smoothed_vorticity_index) * dpsifdx(l, j);
1364 }
1365 } // for (unsigned l=0;l<n_node;l++)
1366 } // End of get_raw_vorticity_deriv
1367
1368 /// Get raw derivative of smoothed vorticity
1370 double& dvorticity_dx,
1371 const unsigned& index) const
1372 {
1373 // Find out how many nodes there are
1374 unsigned n_node = this->nnode();
1375
1376 // Set up memory for the shape functions
1377 Shape psif(n_node);
1378
1379 // Set up memory for the shape function derivatives
1380 DShape dpsifdx(n_node, 2);
1381
1382 // Call the derivatives of the shape and test functions
1383 this->dshape_eulerian(s, psif, dpsifdx);
1384
1385 // Initialise value to zero
1386 dvorticity_dx = 0.0;
1387
1388 // Use a switch statement
1389 switch (index)
1390 {
1391 case 0:
1392 // Loop over the nodes
1393 for (unsigned l = 0; l < n_node; l++)
1394 {
1395 // Calculate the x derivative
1396 dvorticity_dx +=
1397 this->nodal_value(l, Smoothed_vorticity_index) * dpsifdx(l, 0);
1398 }
1399 break;
1400 case 1:
1401 // Loop over the nodes
1402 for (unsigned l = 0; l < n_node; l++)
1403 {
1404 // Calculate the y derivative
1405 dvorticity_dx +=
1406 this->nodal_value(l, Smoothed_vorticity_index) * dpsifdx(l, 1);
1407 }
1408 break;
1409 default:
1410 oomph_info << "Never get here!" << std::endl;
1411 abort();
1412 }
1413 } // End of get_raw_vorticity_deriv
1414
1415 /// Get raw derivative of smoothed derivative vorticity
1417 Vector<double>& dvorticity_dxdy) const
1418 {
1419 // Find out how many nodes there are
1420 unsigned n_node = this->nnode();
1421
1422 // Set up memory for the shape functions
1423 Shape psif(n_node);
1424
1425 // Set up memory for the shape function derivatives
1426 DShape dpsifdx(n_node, 2);
1427
1428 // Call the derivatives of the shape and test functions
1429 this->dshape_eulerian(s, psif, dpsifdx);
1430
1431 // Initialise all entries to zero
1432 dvorticity_dxdy.initialise(0.0);
1433
1434 // Loop over nodes
1435 for (unsigned l = 0; l < n_node; l++)
1436 {
1437 // Loop over derivative directions to obtain xx and xy derivatives
1438 for (unsigned j = 0; j < 2; j++)
1439 {
1440 // Calculate xx and xy derivative
1441 dvorticity_dxdy[j] +=
1442 this->nodal_value(l, Smoothed_vorticity_index + 1) * dpsifdx(l, j);
1443 }
1444
1445 // Calculate the yy derivative
1446 dvorticity_dxdy[2] +=
1447 this->nodal_value(l, Smoothed_vorticity_index + 2) * dpsifdx(l, 1);
1448 } // for (unsigned l=0;l<n_node;l++)
1449 } // End of get_raw_vorticity_second_deriv
1450
1451
1452 /// Get raw derivative of smoothed derivative vorticity
1453 /// [0]: d^2/dx^2, [1]: d^2/dxdy, [2]: d^2/dy^2
1455 double& dvorticity_dxdy,
1456 const unsigned& index) const
1457 {
1458 // Find out how many nodes there are
1459 unsigned n_node = this->nnode();
1460
1461 // Set up memory for the shape functions
1462 Shape psif(n_node);
1463
1464 // Set up memory for the shape function derivatives
1465 DShape dpsifdx(n_node, 2);
1466
1467 // Call the derivatives of the shape and test functions
1468 this->dshape_eulerian(s, psif, dpsifdx);
1469
1470 // Initialise value to zero
1471 dvorticity_dxdy = 0.0;
1472
1473 // Use a switch statement
1474 switch (index)
1475 {
1476 case 0:
1477 // Loop over the nodes
1478 for (unsigned l = 0; l < n_node; l++)
1479 {
1480 // Calculate xx derivative
1481 dvorticity_dxdy +=
1482 this->nodal_value(l, Smoothed_vorticity_index + 1) *
1483 dpsifdx(l, 0);
1484 }
1485 break;
1486 case 1:
1487 // Loop over the nodes
1488 for (unsigned l = 0; l < n_node; l++)
1489 {
1490 // Calculate xy derivative
1491 dvorticity_dxdy +=
1492 this->nodal_value(l, Smoothed_vorticity_index + 1) *
1493 dpsifdx(l, 1);
1494 }
1495 break;
1496 case 2:
1497 // Loop over the nodes
1498 for (unsigned l = 0; l < n_node; l++)
1499 {
1500 // Calculate the yy derivative
1501 dvorticity_dxdy +=
1502 this->nodal_value(l, Smoothed_vorticity_index + 2) *
1503 dpsifdx(l, 1);
1504 }
1505 break;
1506 default:
1507 oomph_info << "Never get here!" << std::endl;
1508 abort();
1509 }
1510 } // End of get_raw_vorticity_second_deriv
1511
1512 /// Get raw derivative of smoothed derivative vorticity
1513 /// [0]: d^3/dx^3, [1]: d^3/dx^2dy, [2]: d^3/dxdy^2, [3]: d^3/dy^3
1515 Vector<double>& dvorticity_dxdxdy) const
1516 {
1517 // Find out how many nodes there are
1518 unsigned n_node = this->nnode();
1519
1520 // Set up memory for the shape functions
1521 Shape psif(n_node);
1522
1523 // Set up memory for the shape function derivatives
1524 DShape dpsifdx(n_node, 2);
1525
1526 // Call the derivatives of the shape and test functions
1527 this->dshape_eulerian(s, psif, dpsifdx);
1528
1529 // Initialise all entries to zero
1530 dvorticity_dxdxdy.initialise(0.0);
1531
1532 // Loop over the nodes
1533 for (unsigned l = 0; l < n_node; l++)
1534 {
1535 // d^3/dx^3 = d/dx \overline{d^2/dx^2}
1536 dvorticity_dxdxdy[0] +=
1537 this->nodal_value(l, Smoothed_vorticity_index + 3) * dpsifdx(l, 0);
1538
1539 // d^3/dx^2dy = d/dx \overline{d^2/dxdy}
1540 dvorticity_dxdxdy[1] +=
1541 this->nodal_value(l, Smoothed_vorticity_index + 4) * dpsifdx(l, 0);
1542
1543 // d^3/dxdy^2 = d/dy \overline{d^2/dxdy}
1544 dvorticity_dxdxdy[2] +=
1545 this->nodal_value(l, Smoothed_vorticity_index + 4) * dpsifdx(l, 1);
1546
1547 // d^3/dy^3 = d/dy \overline{d^2/dy^2}
1548 dvorticity_dxdxdy[3] +=
1549 this->nodal_value(l, Smoothed_vorticity_index + 5) * dpsifdx(l, 1);
1550 }
1551 } // End of get_raw_vorticity_third_deriv
1552
1553
1554 /// Get raw derivative of smoothed derivative vorticity
1555 /// [0]: d^3/dx^3, [1]: d^3/dx^2dy, [2]: d^3/dxdy^2, [3]: d^3/dy^3,
1557 double& dvorticity_dxdxdy,
1558 const unsigned& index) const
1559 {
1560 // Find out how many nodes there are
1561 unsigned n_node = this->nnode();
1562
1563 // Set up memory for the shape functions
1564 Shape psif(n_node);
1565
1566 // Set up memory for the shape function derivatives
1567 DShape dpsifdx(n_node, 2);
1568
1569 // Call the derivatives of the shape and test functions
1570 this->dshape_eulerian(s, psif, dpsifdx);
1571
1572 // Initialise value to zero
1573 dvorticity_dxdxdy = 0.0;
1574
1575 // Use a switch statement
1576 switch (index)
1577 {
1578 case 0:
1579 // Loop over the nodes
1580 for (unsigned l = 0; l < n_node; l++)
1581 {
1582 // d^3/dx^3 = d/dx \overline{d^2/dx^2}
1583 dvorticity_dxdxdy +=
1584 this->nodal_value(l, Smoothed_vorticity_index + 3) *
1585 dpsifdx(l, 0);
1586 }
1587 break;
1588 case 1:
1589 // Loop over the nodes
1590 for (unsigned l = 0; l < n_node; l++)
1591 {
1592 // d^3/dx^2dy = d/dx \overline{d^2/dxdy}
1593 dvorticity_dxdxdy +=
1594 this->nodal_value(l, Smoothed_vorticity_index + 4) *
1595 dpsifdx(l, 0);
1596 }
1597 break;
1598 case 2:
1599 // Loop over the nodes
1600 for (unsigned l = 0; l < n_node; l++)
1601 {
1602 // d^3/dxdy^2 = d/dy \overline{d^2/dxdy}
1603 dvorticity_dxdxdy +=
1604 this->nodal_value(l, Smoothed_vorticity_index + 4) *
1605 dpsifdx(l, 1);
1606 }
1607 break;
1608 case 3:
1609 // Loop over the nodes
1610 for (unsigned l = 0; l < n_node; l++)
1611 {
1612 // d^3/dy^3 = d/dy \overline{d^2/dy^2}
1613 dvorticity_dxdxdy +=
1614 this->nodal_value(l, Smoothed_vorticity_index + 5) *
1615 dpsifdx(l, 1);
1616 }
1617 break;
1618 default:
1619 oomph_info << "Never get here!" << std::endl;
1620 abort();
1621 }
1622 } // End of get_raw_vorticity_third_deriv
1623
1624 /// Compute the element's contribution to the (squared) L2 norm
1625 /// of the difference between exact and smoothed vorticity. The input
1626 /// i corresponds to the i-th dof stored at each node (excluding the
1627 /// velocities and pressure).
1628 double vorticity_error_squared(const unsigned& i)
1629 {
1630#ifdef PARANOID
1631 // Number of derivatives to be recovered
1632 unsigned n_recovered_derivs = (nvorticity_derivatives_to_recover() +
1634
1635 // We cannot calculate this if we're not recovering the data
1636 if ((n_recovered_derivs == 0) || (i >= n_recovered_derivs))
1637 {
1638 // Throw an error
1639 throw OomphLibError("Can't calculate this; not recovering enough data.",
1640 OOMPH_CURRENT_FUNCTION,
1641 OOMPH_EXCEPTION_LOCATION);
1642 }
1643#endif
1644
1645 // Create a container for the synthetic quantities
1646 Vector<Vector<double>> synth_vort_and_derivs =
1648
1649 // Get the appropriate indices associated with the i-th recovered dof
1650 std::pair<unsigned, unsigned> indices = recovered_dof_to_container_id(i);
1651
1652 // Find out how much information we need to calculate
1653 unsigned n_vector = synth_vort_and_derivs.size();
1654
1655 // Initialise the norm squared value
1656 double norm_squared = 0.0;
1657
1658 // Find out how many nodes there are in the element
1659 unsigned n_node = this->nnode();
1660
1661 // Set up memory for the shape functions
1662 Shape psif(n_node);
1663
1664 // Set up memory for the shape function derivatives
1665 DShape dpsifdx(n_node, 2);
1666
1667 // Number of integration points
1668 unsigned n_intpt = this->integral_pt()->nweight();
1669
1670 // Create the vector to hold local coordinates
1671 Vector<double> s(N_dim, 0.0);
1672
1673 // Create the vector to hold the global coordinates
1674 Vector<double> x(N_dim, 0.0);
1675
1676 // Loop over the integration points
1677 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1678 {
1679 // Loop over the coordinates
1680 for (unsigned ii = 0; ii < N_dim; ii++)
1681 {
1682 // Get the local coordinate value
1683 s[ii] = this->integral_pt()->knot(ipt, ii);
1684 }
1685
1686 // Calculate the corresponding global coordinate
1687 this->interpolated_x(s, x);
1688
1689 // Get the integral weight
1690 double w = this->integral_pt()->weight(ipt);
1691
1692 // Call the derivatives of the shape and test functions
1693 double J = this->dshape_eulerian(s, psif, dpsifdx);
1694
1695 // Pre-multiply the weights and the Jacobian
1696 double W = w * J;
1697
1698 // Initialise the smoothed vorticity value
1699 double smoothed_vort = 0.0;
1700
1701 // Smoothed vorticity
1702 for (unsigned l = 0; l < n_node; l++)
1703 {
1704 // Update the smoothed vorticity value
1705 smoothed_vort +=
1706 this->nodal_value(l, Smoothed_vorticity_index + i) * psif[l];
1707 }
1708
1709 // Initialise the entries of the storage for the vorticity and
1710 // derivatives
1711 for (unsigned jj = 0; jj < n_vector; jj++)
1712 {
1713 // Initialise the entries to zero
1714 synth_vort_and_derivs[jj].initialise(0.0);
1715 }
1716
1717 // If pointer isn't null
1718 if (0 != Exact_vorticity_fct_pt)
1719 {
1720 // Use the function pointer to calculate the appropriate quantities
1721 Exact_vorticity_fct_pt(x, synth_vort_and_derivs);
1722 }
1723
1724 // Initialise the synthetic quantity
1725 double synth_quantity = 0.0;
1726
1727 // Calculate the synthetic quantity
1728 synth_quantity = (synth_vort_and_derivs[indices.first])[indices.second];
1729
1730 // Add the squared difference
1731 norm_squared += pow(smoothed_vort - synth_quantity, 2) * W;
1732 } // for (unsigned ipt=0;ipt<n_intpt;ipt++)
1733
1734 // Return the norm squared value
1735 return norm_squared;
1736 } // End of vorticity_error_squared
1737
1738
1739 /// Compute smoothed vorticity and its derivatives
1741 Vector<Vector<double>>& vort_and_derivs) const
1742 {
1743 // Get the number of nodes in this element
1744 unsigned n_node = this->nnode();
1745
1746 // Shape functions
1747 Shape psif(n_node);
1748
1749 // Get the shape function value at the given local coordinate
1750 this->shape(s, psif);
1751
1752 // Find out how much information we need to calculate
1753 unsigned n_vector = vort_and_derivs.size();
1754
1755 // Initialise a counter (holds the dof number)
1756 unsigned i_dof = 0;
1757
1758 // Loop over the vectors
1759 for (unsigned i = 0; i < n_vector; i++)
1760 {
1761 // Initialise entries to zero
1762 vort_and_derivs[i].initialise(0.0);
1763
1764 // Find out how many entries there are in the i-th vector
1765 unsigned num_entries = vort_and_derivs[i].size();
1766
1767 // Loop over the entries
1768 for (unsigned j = 0; j < num_entries; j++)
1769 {
1770 // Loop over the nodes
1771 for (unsigned l = 0; l < n_node; l++)
1772 {
1773 // Get the contribution to the smoothed derivative from the l-th
1774 // node
1775 (vort_and_derivs[i])[j] +=
1776 this->nodal_value(l, Smoothed_vorticity_index + i_dof) * psif[l];
1777 }
1778
1779 // Increment the counter
1780 i_dof++;
1781 } // for (unsigned j=0;j<num_entries;j++)
1782 } // for (unsigned i=0;i<n_vector;i++)
1783 } // End of vorticity_and_its_derivs
1784
1785 private:
1786 /// Number of dimensions in the element
1787 unsigned N_dim;
1788
1789 /// Index of smoothed vorticity -- followed by derivatives;
1790 /// in 2D this has value 3
1792
1793 /// The current maximum order of vorticity derivatives that can be
1794 /// recovered. Currently, we can recover up to the third derivative:
1795 /// omega,d/dx,d/dy,
1796 /// d^2/dx^2,d^2/dxdy,d^2/dy^2,
1797 /// d^3/dx^3,d^3/dx^2dy,d^3/dxdy^2,d^3/dy^3
1799
1800 /// The current maximum order of velocity derivatives that can be
1801 /// recovered. Currently, we can recover the first derivatives:
1802 /// du/dx,du/dy,dv/dx,dv/dy
1804
1805 /// Number of values per field; how many of the following do we want:
1806 /// u,v,p,omega,d/dx,d/dy,
1807 /// d^2/dx^2,d^2/dxdy,d^2/dy^2,
1808 /// d^3/dx^3,d^3/dx^2dy,d^3/dxdy^2,d^3/dy^3,
1809 /// du/dx,du/dy,dv/dx,dv/dy
1811
1812 /// Maximum number of derivatives to retain in the vorticity
1813 /// recovery. Note, the value -1 means we ONLY output u,v[,w],p.
1815
1816 /// Maximum number of derivatives to retain in the velocity
1817 /// recovery. Note, the value 0 means we don't calculate the derivatives
1818 /// of the velocity
1820
1821 /// Pointer to function that specifies exact vorticity and
1822 /// derivs (for validation).
1824 };
1825
1826
1827 /// ///////////////////////////////////////////////////////////////////
1828 /// ///////////////////////////////////////////////////////////////////
1829 /// ///////////////////////////////////////////////////////////////////
1830
1831
1832 //========================================================
1833 /// Smoother for vorticity in 2D
1834 //========================================================
1835 template<class ELEMENT>
1837 {
1838 public:
1839 /// Constructor: Set order of recovery shape functions
1842 {
1843 }
1844
1845 /// Broken copy constructor
1847
1848 /// Broken assignment operator
1849 void operator=(const VorticitySmoother&) = delete;
1850
1851 /// Empty virtual destructor
1853
1854 /// Access function for order of recovery polynomials
1855 unsigned& recovery_order()
1856 {
1857 // Return the order of recovery
1858 return Recovery_order;
1859 }
1860
1861 /// Recovery shape functions as functions of the global, Eulerian
1862 /// coordinate x of dimension dim. The recovery shape functions are complete
1863 /// polynomials of the order specified by Recovery_order.
1865 {
1866 // Create an ostringstream object to create a string
1867 std::ostringstream error_stream;
1868
1869 // Find order of recovery shape functions
1870 switch (recovery_order())
1871 {
1872 case 1:
1873 // Complete linear polynomial in 2D:
1874 psi_r[0] = 1.0;
1875 psi_r[1] = x[0];
1876 psi_r[2] = x[1];
1877 break;
1878
1879 case 2:
1880 // Complete quadratic polynomial in 2D:
1881 psi_r[0] = 1.0;
1882 psi_r[1] = x[0];
1883 psi_r[2] = x[1];
1884 psi_r[3] = x[0] * x[0];
1885 psi_r[4] = x[0] * x[1];
1886 psi_r[5] = x[1] * x[1];
1887 break;
1888
1889 case 3:
1890 // Complete cubic polynomial in 2D:
1891 psi_r[0] = 1.0;
1892 psi_r[1] = x[0];
1893 psi_r[2] = x[1];
1894 psi_r[3] = x[0] * x[0];
1895 psi_r[4] = x[0] * x[1];
1896 psi_r[5] = x[1] * x[1];
1897 psi_r[6] = x[0] * x[0] * x[0];
1898 psi_r[7] = x[0] * x[0] * x[1];
1899 psi_r[8] = x[0] * x[1] * x[1];
1900 psi_r[9] = x[1] * x[1] * x[1];
1901 break;
1902
1903 default:
1904 // Create an error message for this case
1905 error_stream << "Recovery shape functions for recovery order "
1906 << recovery_order()
1907 << " haven't yet been implemented for 2D" << std::endl;
1908
1909 // Throw an error
1910 throw OomphLibError(error_stream.str(),
1911 OOMPH_CURRENT_FUNCTION,
1912 OOMPH_EXCEPTION_LOCATION);
1913 }
1914 } // End of shape_rec
1915
1916
1917 /// Integation scheme associated with the recovery shape functions
1918 /// must be of sufficiently high order to integrate the mass matrix
1919 /// associated with the recovery shape functions. The argument is the
1920 /// dimension of the elements.
1921 /// The integration is performed locally over the elements, so the
1922 /// integration scheme does depend on the geometry of the element.
1923 /// The type of element is specified by the boolean which is
1924 /// true if elements in the patch are QElements and false if they are
1925 /// TElements (will need change if we ever have other element types)
1926 Integral* integral_rec(const bool& is_q_mesh)
1927 {
1928 // Create an ostringstream object to create a string
1929 std::ostringstream error_stream;
1930
1931 //----
1932 // 2D:
1933 //----
1934 /// Find order of recovery shape functions
1935 switch (recovery_order())
1936 {
1937 case 1:
1938 // Complete linear polynomial in 2D:
1939 if (is_q_mesh)
1940 {
1941 // Return the appropriate Gauss integration scheme
1942 return (new Gauss<2, 2>);
1943 }
1944 else
1945 {
1946 // Return the appropriate Gauss integration scheme
1947 return (new TGauss<2, 2>);
1948 }
1949 break;
1950
1951 case 2:
1952 // Complete quadratic polynomial in 2D:
1953 if (is_q_mesh)
1954 {
1955 // Return the appropriate Gauss integration scheme
1956 return (new Gauss<2, 3>);
1957 }
1958 else
1959 {
1960 // Return the appropriate Gauss integration scheme
1961 return (new TGauss<2, 3>);
1962 }
1963 break;
1964
1965 case 3:
1966 // Complete cubic polynomial in 2D:
1967 if (is_q_mesh)
1968 {
1969 // Return the appropriate Gauss integration scheme
1970 return (new Gauss<2, 4>);
1971 }
1972 else
1973 {
1974 // Return the appropriate Gauss integration scheme
1975 return (new TGauss<2, 4>);
1976 }
1977 break;
1978
1979 default:
1980 // Create an error messaage
1981 error_stream << "Recovery shape functions for recovery order "
1982 << recovery_order()
1983 << " haven't yet been implemented for 2D" << std::endl;
1984
1985 // Throw an error
1986 throw OomphLibError(error_stream.str(),
1987 OOMPH_CURRENT_FUNCTION,
1988 OOMPH_EXCEPTION_LOCATION);
1989 }
1990
1991 // Dummy return (never get here)
1992 return 0;
1993 } // End of integral_rec
1994
1995 /// Setup patches: For each vertex node pointed to by nod_pt,
1996 /// adjacent_elements_pt[nod_pt] contains the pointer to the vector that
1997 /// contains the pointers to the elements that the node is part of.
1998 /// Also returns a Vector of vertex nodes for use in get_element_errors.
1999 void setup_patches(Mesh*& mesh_pt,
2000 std::map<Node*, Vector<ELEMENT*>*>& adjacent_elements_pt,
2001 Vector<Node*>& vertex_node_pt)
2002 {
2003 // Clear: hierher should we do this in Z2 as well?
2004 adjacent_elements_pt.clear();
2005
2006 // Auxiliary map that contains element-adjacency for ALL nodes
2007 std::map<Node*, Vector<ELEMENT*>*> aux_adjacent_elements_pt;
2008
2009#ifdef PARANOID
2010 // Check if all elements request the same recovery order
2011 unsigned ndisagree = 0;
2012#endif
2013
2014 // Loop over all elements to setup adjacency for all nodes.
2015 // Need to do this because midside nodes can be corner nodes for
2016 // adjacent smaller elements! Admittedly, the inclusion of interior
2017 // nodes is wasteful...
2018 unsigned nelem = mesh_pt->nelement();
2019 for (unsigned e = 0; e < nelem; e++)
2020 {
2021 ELEMENT* el_pt = dynamic_cast<ELEMENT*>(mesh_pt->element_pt(e));
2022
2023#ifdef PARANOID
2024 // Check if all elements request the same recovery order
2025 if (el_pt->nrecovery_order() != Recovery_order)
2026 {
2027 ndisagree++;
2028 }
2029#endif
2030
2031 // Loop all nodes in element
2032 unsigned nnod = el_pt->nnode();
2033 for (unsigned n = 0; n < nnod; n++)
2034 {
2035 // Make a pointer to the n-th node
2036 Node* nod_pt = el_pt->node_pt(n);
2037
2038 // Has this node been considered before?
2039 if (aux_adjacent_elements_pt[nod_pt] == 0)
2040 {
2041 // Create Vector of pointers to its adjacent elements
2042 aux_adjacent_elements_pt[nod_pt] = new Vector<ELEMENT*>;
2043 }
2044
2045 // Add pointer to adjacent element
2046 (*aux_adjacent_elements_pt[nod_pt]).push_back(el_pt);
2047 }
2048 } // end element loop
2049
2050#ifdef PARANOID
2051 // Check if all elements request the same recovery order
2052 if (ndisagree != 0)
2053 {
2055 << "\n\n======================================================\n"
2056 << "WARNING:\n"
2057 << ndisagree << " out of " << mesh_pt->nelement() << " elements"
2058 << "\nhave different preferences for the order of the recovery"
2059 << "\nshape functions. We are using: Recovery_order="
2060 << Recovery_order << std::endl;
2062 << "======================================================\n\n";
2063 }
2064#endif
2065
2066 // Loop over all elements, extract adjacency for corner nodes only
2067 nelem = mesh_pt->nelement();
2068 for (unsigned e = 0; e < nelem; e++)
2069 {
2070 ELEMENT* el_pt = dynamic_cast<ELEMENT*>(mesh_pt->element_pt(e));
2071
2072 // Loop over corner nodes
2073 unsigned n_node = el_pt->nvertex_node();
2074 for (unsigned n = 0; n < n_node; n++)
2075 {
2076 Node* nod_pt = el_pt->vertex_node_pt(n);
2077
2078 // Has this node been considered before?
2079 if (adjacent_elements_pt[nod_pt] == 0)
2080 {
2081 // Add the node pointer to the vertex node container
2082 vertex_node_pt.push_back(nod_pt);
2083
2084 // Create Vector of pointers to its adjacent elements
2085 adjacent_elements_pt[nod_pt] = new Vector<ELEMENT*>;
2086
2087 // Copy across:
2088 unsigned nel = (*aux_adjacent_elements_pt[nod_pt]).size();
2089 for (unsigned e = 0; e < nel; e++)
2090 {
2091 (*adjacent_elements_pt[nod_pt])
2092 .push_back((*aux_adjacent_elements_pt[nod_pt])[e]);
2093 }
2094 }
2095 }
2096 } // End of loop over elements
2097
2098 // Cleanup
2099 for (typename std::map<Node*, Vector<ELEMENT*>*>::iterator it =
2100 aux_adjacent_elements_pt.begin();
2101 it != aux_adjacent_elements_pt.end();
2102 it++)
2103 {
2104 delete it->second;
2105 }
2106 } // End of setup_patches
2107
2108 /// Given the vector of elements that make up a patch, compute
2109 /// the vector of recovered vorticity coefficients and return a pointer
2110 /// to it. n_deriv indicates which derivative of the vorticity is
2111 /// supposed to be smoothed: 0: zeroth (i.e. the vorticity itself)
2112 /// 1: d/dx; 2: d/dy; 3: d^2/dx^2; 4: d^2/dxdy 5: d^2/dy^2
2113 /// 6: d^3/dx^3, 7: d^3/dx^2dy, 8: d^3/dxdy^2, 9: d^3/dy^3,
2114 /// 10: du/dx, 11: du/dy, 12: dv/dx, 13: dv/dy
2116 const Vector<ELEMENT*>& patch_el_pt,
2117 const unsigned& num_recovery_terms,
2118 Vector<double>*& recovered_vorticity_coefficient_pt,
2119 unsigned& n_deriv)
2120 {
2121 // Find the number of elements in the patch
2122 unsigned nelem = patch_el_pt.size();
2123
2124 // Get a pointer to any element
2125 ELEMENT* el_pt = patch_el_pt[0];
2126
2127#ifdef PARANOID
2128 // If there's at least one element
2129 if (nelem > 0)
2130 {
2131 // Get the number of vorticity derivatives to recover
2132 int n_vort_derivs = el_pt->get_maximum_order_of_vorticity_derivative();
2133
2134 // Get the number of vorticity derivatives to recover
2135 int n_veloc_derivs = el_pt->get_maximum_order_of_velocity_derivative();
2136
2137 // If we're not recovering anything, we shouldn't be here
2138 if (n_vort_derivs + n_veloc_derivs == -1)
2139 {
2140 // Create an ostringstream object to create an error message
2141 std::ostringstream error_stream;
2142
2143 // Create the error message
2144 error_stream << "Not recovering anything. Change the maximum number "
2145 << "of derivatives to recover.";
2146
2147 // Throw an error
2148 throw OomphLibError(error_stream.str(),
2149 OOMPH_CURRENT_FUNCTION,
2150 OOMPH_EXCEPTION_LOCATION);
2151 }
2152 } // if (nelem>0)
2153#endif
2154
2155 // Find the container indices associated with n_deriv
2156 std::pair<unsigned, unsigned> container_id =
2157 el_pt->vorticity_dof_to_container_id(n_deriv);
2158
2159 // Maximum vorticity derivative order we can recover
2160 unsigned max_recoverable_vort_order =
2161 el_pt->get_maximum_order_of_recoverable_vorticity_derivative();
2162
2163 // Maximum velocity derivative order we can recover
2164 unsigned max_recoverable_veloc_order =
2165 el_pt->get_maximum_order_of_recoverable_velocity_derivative();
2166
2167 // Make a counter
2168 unsigned counter = 0;
2169
2170 // Calculate the case value (initialise to -1 so we know if it's set
2171 // later)
2172 int case_value = -1;
2173
2174 // Loop over the derivatives
2175 for (unsigned i = 0; i < max_recoverable_vort_order + 1; i++)
2176 {
2177 // Increment by the number of partial derivatives of order i
2178 counter += el_pt->npartial_derivative(i);
2179
2180 // If we've exceeded the value of n_deriv then we know which vorticity
2181 // derivative to recover
2182 if (n_deriv < counter)
2183 {
2184 // We need to recover the i-th order of derivative of the vorticity
2185 case_value = i;
2186
2187 // We're done here
2188 break;
2189 }
2190 } // for (unsigned i=0;i<max_recoverable_order+1;i++)
2191
2192 // If we haven't set the case value yet then we must be recovering a
2193 // velocity derivative
2194 if (case_value == -1)
2195 {
2196 // Loop over the velocity order
2197 for (unsigned i = 1; i < max_recoverable_veloc_order + 1; i++)
2198 {
2199 // Increment by the number of velocity partial derivatives of order i
2200 counter += 2 * el_pt->npartial_derivative(i);
2201
2202 // If we've exceeded the value of n_deriv then we know which vorticity
2203 // derivative to recover
2204 if (n_deriv < counter)
2205 {
2206 // We need to recover the i-th order of derivative of the vorticity
2207 case_value = max_recoverable_vort_order + i;
2208
2209 // We're done here
2210 break;
2211 }
2212 } // for (unsigned i=1;i<max_recoverable_veloc_order+1;i++)
2213 } // if (case_value==-1)
2214
2215#ifdef PARANOID
2216 // Sanity check: if the case value hasn't been set then something's wrong
2217 if (case_value == -1)
2218 {
2219 // Create a ostringstream object to create an error message
2220 std::ostringstream error_message_stream;
2221
2222 // Create an error message
2223 error_message_stream
2224 << "Case order has not been set. Something's wrong!";
2225
2226 // Throw an error
2227 throw OomphLibError(error_message_stream.str(),
2228 OOMPH_CURRENT_FUNCTION,
2229 OOMPH_EXCEPTION_LOCATION);
2230 }
2231#endif
2232
2233 // Create space for the recovered quantity
2234 double recovered_quantity = 0.0;
2235
2236 // Create/initialise matrix for linear system
2237 DenseDoubleMatrix recovery_mat(
2238 num_recovery_terms, num_recovery_terms, 0.0);
2239
2240 // Create/initialise vector for RHS
2241 Vector<double> rhs(num_recovery_terms, 0.0);
2242
2243 // Create a new integration scheme based on the recovery order
2244 // in the elements. Need to find the type of the element, default
2245 // is to assume a quad
2246 bool is_q_mesh = true;
2247
2248 // If we can dynamic cast to the TElementBase, then it's a triangle/tet
2249 // Note that I'm assuming that all elements are of the same geometry, but
2250 // if they weren't we could adapt...
2251 if (dynamic_cast<TElementBase*>(patch_el_pt[0]))
2252 {
2253 // We're dealing with a triangle-based mesh so change the bool value
2254 is_q_mesh = false;
2255 }
2256
2257 // Get a pointer to the appropriate integration type
2258 Integral* const integ_pt = this->integral_rec(is_q_mesh);
2259
2260 // Loop over all elements in patch to assemble linear system
2261 for (unsigned e = 0; e < nelem; e++)
2262 {
2263 // Get pointer to element
2264 ELEMENT* const el_pt = patch_el_pt[e];
2265
2266 // Create storage for the recovery shape function values
2267 Vector<double> psi_r(num_recovery_terms);
2268
2269 // Create vector to hold local coordinates
2270 Vector<double> s(2, 0.0);
2271
2272 // Find the number of integration points
2273 unsigned n_intpt = integ_pt->nweight();
2274
2275 // Loop over the integration points
2276 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
2277 {
2278 // Assign values of s, the local coordinate
2279 for (unsigned i = 0; i < 2; i++)
2280 {
2281 // Get the i-th entry of the local coordinate
2282 s[i] = integ_pt->knot(ipt, i);
2283 }
2284
2285 // Get the integral weight
2286 double w = integ_pt->weight(ipt);
2287
2288 // Jacobian of mapping
2289 double J = el_pt->J_eulerian(s);
2290
2291 // Allocate space for the global coordinates
2292 Vector<double> x(2, 0.0);
2293
2294 // Interpolate the global (Eulerian) coordinate
2295 el_pt->interpolated_x(s, x);
2296
2297 // Premultiply the weights and the Jacobian and the geometric
2298 // Jacobian weight (used in axisymmetric and spherical coordinate
2299 // systems) -- hierher really fct of x? probably yes, actually).
2300 double W = w * J * (el_pt->geometric_jacobian(x));
2301
2302 // Recovery shape functions at global (Eulerian) coordinate
2303 shape_rec(x, psi_r);
2304
2305 // Use a switch statement to decide on which function to call
2306 switch (case_value)
2307 {
2308 case 0:
2309 {
2310 Vector<double> vorticity(1, 0.0);
2311 el_pt->get_vorticity(s, vorticity);
2312 recovered_quantity = vorticity[0];
2313 break;
2314 }
2315 case 1:
2316 {
2317 el_pt->get_raw_vorticity_deriv(
2318 s, recovered_quantity, container_id.second);
2319 break;
2320 }
2321 case 2:
2322 {
2323 el_pt->get_raw_vorticity_second_deriv(
2324 s, recovered_quantity, container_id.second);
2325 break;
2326 }
2327 case 3:
2328 {
2329 el_pt->get_raw_vorticity_third_deriv(
2330 s, recovered_quantity, container_id.second);
2331 break;
2332 }
2333 case 4:
2334 {
2335 el_pt->get_raw_velocity_deriv(
2336 s, recovered_quantity, container_id.second);
2337 break;
2338 }
2339 default:
2340 {
2341 oomph_info << "Never get here." << std::endl;
2342 abort();
2343 }
2344 }
2345
2346 // Add elemental RHSs and recovery matrix to global versions
2347 //----------------------------------------------------------
2348 // RHS: Loop over the nodes for the test functions
2349 for (unsigned l = 0; l < num_recovery_terms; l++)
2350 {
2351 // Update the RHS entry ()
2352 rhs[l] += recovered_quantity * psi_r[l] * W;
2353 }
2354
2355 // Loop over the nodes for the test functions
2356 for (unsigned l = 0; l < num_recovery_terms; l++)
2357 {
2358 // Loop over the nodes for the variables
2359 for (unsigned l2 = 0; l2 < num_recovery_terms; l2++)
2360 {
2361 // Add contribution to recovery matrix
2362 recovery_mat(l, l2) += psi_r[l] * psi_r[l2] * W;
2363 }
2364 } // for (unsigned l=0;l<num_recovery_terms;l++)
2365 } // for (unsigned ipt=0;ipt<n_intpt;ipt++)
2366 } // End of loop over elements that make up patch.
2367
2368 // Delete the integration scheme
2369 delete integ_pt;
2370
2371 // Linear system is now assembled: Solve recovery system
2372 //------------------------------------------------------
2373 // LU decompose the recovery matrix
2374 recovery_mat.ludecompose();
2375
2376 // Back-substitute
2377 recovery_mat.lubksub(rhs);
2378
2379 // Now create a matrix to store the vorticity recovery coefficients.
2380 // Pointer to this matrix will be returned.
2381 recovered_vorticity_coefficient_pt =
2382 new Vector<double>(num_recovery_terms);
2383
2384 // Loop over the number of recovered terms
2385 for (unsigned icoeff = 0; icoeff < num_recovery_terms; icoeff++)
2386 {
2387 // Copy the RHS value over
2388 (*recovered_vorticity_coefficient_pt)[icoeff] = rhs[icoeff];
2389 }
2390 } // End of get_recovered_vorticity_in_patch
2391
2392 /// Get the recovery order
2393 unsigned nrecovery_order() const
2394 {
2395 // Use a switch statement
2396 switch (Recovery_order)
2397 {
2398 case 1:
2399 // Linear recovery shape functions
2400 //--------------------------------
2401 return 3; // 1, x, y
2402 break;
2403
2404 case 2:
2405 // Quadratic recovery shape functions
2406 //-----------------------------------
2407 return 6; // 1, x, y, x^2, xy, y^2
2408 break;
2409
2410 case 3:
2411 // Cubic recovery shape functions
2412 //--------------------------------
2413 return 10; // 1, x, y, x^2, xy, y^2, x^3, x^2 y, x y^2, y^3
2414 break;
2415
2416 default:
2417 // Any other recovery order?
2418 //--------------------------
2419 // Use an ostringstream object to create an error message
2420 std::ostringstream error_stream;
2421
2422 // Create an error message
2423 error_stream << "Wrong Recovery_order " << Recovery_order
2424 << std::endl;
2425
2426 // Throw an error
2427 throw OomphLibError(error_stream.str(),
2428 OOMPH_CURRENT_FUNCTION,
2429 OOMPH_EXCEPTION_LOCATION);
2430 }
2431 } // End of nrecovery_order
2432
2433 /// Recover vorticity from patches
2435 {
2436 // Create a DocInfo object (used as a dummy argument)
2437 DocInfo doc_info;
2438
2439 // Disable any documentation
2440 doc_info.disable_doc();
2441
2442 // Recover the vorticity
2443 recover_vorticity(mesh_pt, doc_info);
2444 }
2445
2446 /// Recover vorticity from patches -- output intermediate steps
2447 /// to directory specified by DocInfo object
2448 void recover_vorticity(Mesh* mesh_pt, DocInfo& doc_info)
2449 {
2450 // Start the timer
2451 double t_start = TimingHelpers::timer();
2452
2453 // Allocate space for the local coordinates
2454 Vector<double> s(2, 0.0);
2455
2456 // Allocate space for the global coordinates
2457 Vector<double> x(2, 0.0);
2458
2459 // Make patches
2460 //-------------
2461 // Allocate space for the mapping from nodes to elements
2462 std::map<Node*, Vector<ELEMENT*>*> adjacent_elements_pt;
2463
2464 // Allocate space for the vertex nodes
2465 Vector<Node*> vertex_node_pt;
2466
2467 // Set up the patches information
2468 setup_patches(mesh_pt, adjacent_elements_pt, vertex_node_pt);
2469
2470 // Grab any element (this shouldn't be a null pointer)
2471 ELEMENT* const el_pt = dynamic_cast<ELEMENT*>(mesh_pt->element_pt(0));
2472
2473 // Get the index of the vorticity
2474 unsigned smoothed_vorticity_index = el_pt->smoothed_vorticity_index();
2475
2476 // Maximum order of vorticity derivative (that can be recovered)
2477 unsigned max_vort_order =
2478 el_pt->get_maximum_order_of_recoverable_vorticity_derivative();
2479
2480 // Maximum order of velocity derivative (that can be recovered)
2481 unsigned max_veloc_order =
2482 el_pt->get_maximum_order_of_recoverable_velocity_derivative();
2483
2484 // Maximum number of recoverable vorticity terms
2485 unsigned max_vort_recov = 0;
2486
2487 // Maximum number of recoverable velocity terms
2488 unsigned max_veloc_recov = 0;
2489
2490 // Loop over the entries of the vector
2491 for (unsigned i = 0; i < max_vort_order + 1; i++)
2492 {
2493 // Get the number of partial derivatives of the vorticity
2494 max_vort_recov += el_pt->npartial_derivative(i);
2495 }
2496
2497 // Loop over the entries of the vector
2498 for (unsigned i = 1; i < max_veloc_order + 1; i++)
2499 {
2500 // Get the number of partial derivatives of the vorticity
2501 max_veloc_recov += 2 * el_pt->npartial_derivative(i);
2502 }
2503
2504 // Number of recovered vorticity derivatives
2505 unsigned n_recovered_vort_derivs =
2506 el_pt->nvorticity_derivatives_to_recover();
2507
2508 // Number of recovered velocity derivatives
2509 unsigned n_recovered_veloc_derivs =
2510 el_pt->nvelocity_derivatives_to_recover();
2511
2512 // Determine number of coefficients for expansion of recovered vorticity.
2513 // Use complete polynomial of given order for recovery
2514 unsigned num_recovery_terms = nrecovery_order();
2515
2516 // Counter for averaging of recovered vorticity and its derivatives
2517 std::map<Node*, unsigned> count;
2518
2519 // Counter for which nodal value we're assigning
2520 unsigned nodal_dof = 0;
2521
2522 // Loop over derivatives
2523 for (unsigned deriv = 0; deriv < max_vort_recov + max_veloc_recov;
2524 deriv++)
2525 {
2526 // If we're not recovering this vorticity derivative. Note, we cast
2527 // to an int because n_recovered_vort_derivs can be zero (so subtracting
2528 // any positive integer can cause trouble)
2529 if ((int(deriv) > int(n_recovered_vort_derivs - 1)) &&
2530 (deriv < max_vort_recov))
2531 {
2532 // We're done here
2533 continue;
2534 }
2535 // If we're not recovering any of the velocity derivatives and we're
2536 // finished with the vorticity derivatives
2537 else if ((n_recovered_veloc_derivs == 0) && (deriv >= max_vort_recov))
2538 {
2539 // We're done here
2540 continue;
2541 }
2542
2543 // Storage for accumulated nodal vorticity (used to compute nodal
2544 // averages)
2545 std::map<Node*, double> averaged_recovered_vort;
2546
2547 // Calculation of vorticity
2548 //-------------------------
2549 // Do patch recovery
2550 // unsigned counter=0;
2551 for (typename std::map<Node*, Vector<ELEMENT*>*>::iterator it =
2552 adjacent_elements_pt.begin();
2553 it != adjacent_elements_pt.end();
2554 it++)
2555 {
2556 // Pointer to the recovered vorticity coefficients
2557 Vector<double>* recovered_vorticity_coefficient_pt;
2558
2559 // Setup smoothed vorticity field for patches
2561 num_recovery_terms,
2562 recovered_vorticity_coefficient_pt,
2563 deriv);
2564
2565 // Now get the nodal average of the recovered vorticity (nodes are
2566 // generally part of multiple patches):
2567
2568 // Get the number of elements in adjacent_elements_pt
2569 unsigned nelem = (*(it->second)).size();
2570
2571 // Loop over all elements to get recovered vorticity
2572 for (unsigned e = 0; e < nelem; e++)
2573 {
2574 // Get pointer to element
2575 ELEMENT* const el_pt = (*(it->second))[e];
2576
2577 // Get the number of nodes by element
2578 unsigned nnode_el = el_pt->nnode();
2579
2580 // Loop over the nodes in the element
2581 for (unsigned j = 0; j < nnode_el; j++)
2582 {
2583 // Get a pointer to the j-th node in this element
2584 Node* nod_pt = el_pt->node_pt(j);
2585
2586 // Get the local coordinates of the node
2587 el_pt->local_coordinate_of_node(j, s);
2588
2589 // Interpolate the global (Eulerian) coordinate
2590 el_pt->interpolated_x(s, x);
2591
2592 // Recovery shape functions at global (Eulerian) coordinate
2593 Vector<double> psi_r(num_recovery_terms);
2594
2595 // Recover the shape function values at the position x
2596 shape_rec(x, psi_r);
2597
2598 // Initialise the value of the recovered quantity
2599 double recovered_vort = 0.0;
2600
2601 // Loop over the recovery terms
2602 for (unsigned i = 0; i < num_recovery_terms; i++)
2603 {
2604 // Assemble recovered vorticity
2605 recovered_vort +=
2606 (*recovered_vorticity_coefficient_pt)[i] * psi_r[i];
2607 }
2608
2609 // Keep adding
2610 averaged_recovered_vort[nod_pt] += recovered_vort;
2611
2612 // Increment the counter
2613 count[nod_pt]++;
2614 } // for (unsigned j=0;j<nnode_el;j++)
2615 } // for (unsigned e=0;e<nelem;e++)
2616
2617 // Delete the recovered coefficient data
2618 delete recovered_vorticity_coefficient_pt;
2619
2620 // Make it a null pointer
2621 recovered_vorticity_coefficient_pt = 0;
2622 } // for (typename std::map<Node*,Vector<ELEMENT*>*>::iterator it=...
2623
2624 // Find out how many nodes there are in the mesh
2625 unsigned nnod = mesh_pt->nnode();
2626
2627 // Loop over all nodes to actually work out the average
2628 for (unsigned j = 0; j < nnod; j++)
2629 {
2630 // Make a pointer to the j-th node
2631 Node* nod_pt = mesh_pt->node_pt(j);
2632
2633 // Calculate the values of the smoothed vorticity
2634 averaged_recovered_vort[nod_pt] /= count[nod_pt];
2635
2636 // Assign smoothed vorticity to nodal values
2637 nod_pt->set_value(smoothed_vorticity_index + nodal_dof,
2638 averaged_recovered_vort[nod_pt]);
2639 }
2640
2641 // We're done with this dof so increment the counter
2642 nodal_dof++;
2643
2644 // Start again
2645 count.clear();
2646 } // for (unsigned deriv=0;deriv<max_vort_recov+max_veloc_recov;deriv++)
2647
2648 // Cleanup
2649 for (typename std::map<Node*, Vector<ELEMENT*>*>::iterator it =
2650 adjacent_elements_pt.begin();
2651 it != adjacent_elements_pt.end();
2652 it++)
2653 {
2654 // Delete the vector of element pointers
2655 delete it->second;
2656 }
2657
2658 // Inform the user
2659 oomph_info << "Time for vorticity recovery [sec]: "
2660 << TimingHelpers::timer() - t_start << std::endl;
2661 } // End of recover_vorticity
2662
2663 private:
2664 /// Order of recovery polynomials
2666 };
2667
2668} // namespace oomph
2669
2670#endif
e
Definition: cfortran.h:571
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:385
void set_value(const unsigned &i, const double &value_)
Set the i-th stored data value to specified value. The only reason that we require an explicit set fu...
Definition: nodes.h:271
Class of matrices containing doubles, and stored as a DenseMatrix<double>, but with solving functiona...
Definition: matrices.h:1271
virtual void lubksub(DoubleVector &rhs)
LU backsubstitution.
Definition: matrices.cc:202
virtual void ludecompose()
LU decomposition using DenseLU (default linea solver)
Definition: matrices.cc:192
Information for documentation of results: Directory and file number to enable output in the form RESL...
void disable_doc()
Disable documentation.
2D Gaussian integration class. 2x2 integration points. This integration scheme can integrate up to th...
Definition: integral.h:297
2D Gaussian integration class. 3x3 integration points. This integration scheme can integrate up to fi...
Definition: integral.h:343
2D Gaussian integration class. 4x4 integration points. This integration scheme can integrate up to se...
Definition: integral.h:389
Generic class for numerical integration schemes:
Definition: integral.h:49
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
A general mesh class.
Definition: mesh.h:67
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
Definition: mesh.h:436
unsigned long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:596
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:448
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:590
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
An OomphLibError object which should be thrown when an run-time error is encountered....
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
Definition: Telements.h:1130
2D Gaussian integration class for linear triangles. Three integration points. This integration scheme...
Definition: integral.h:820
2D Gaussian integration class for quadratic triangles. Seven integration points. This integration sch...
Definition: integral.h:867
2D Gaussian integration class for cubic triangles. Thirteen integration points. This integration sche...
Definition: integral.h:914
void initialise(const _Tp &__value)
Iterate over all values and set to the desired value.
Definition: Vector.h:167
Class to indicate which derivatives of the vorticity/ velocity we want to recover....
int maximum_order_of_velocity_derivative() const
The maximum order of derivatives calculated in the velocity recovery.
unsigned Number_of_values_per_field
Number of values per field; how many of the following do we want: u,v,p,omega,d/dx,...
int Maximum_order_of_vorticity_derivative
Maximum number of derivatives to retain in the vorticity recovery. Note, the value -1 means we ONLY o...
void set_maximum_order_of_velocity_derivative(const int &max_deriv)
The maximum order of derivatives calculated in the velocity recovery.
void calculate_number_of_values_per_field()
Calculates the number of values per field given the number of vorticity and velocity derivatives to r...
int Maximum_order_of_velocity_derivative
Maximum number of derivatives to retain in the velocity recovery. Note, the value 0 means we don't ca...
unsigned ncont_interpolated_values() const
Number of continuously interpolated values:
unsigned npartial_derivative(const unsigned &n) const
Helper function that determines the number of n-th order partial derivatives in d-dimensions....
void set_maximum_order_of_vorticity_derivative(const int &max_deriv)
The maximum order of derivatives calculated in the vorticity recovery.
int maximum_order_of_vorticity_derivative() const
The maximum order of derivatives calculated in the vorticity recovery.
Overloaded element that allows projection of vorticity.
int get_maximum_order_of_velocity_derivative() const
The maximum order of derivatives calculated in the velocity recovery. Note, this value can only be se...
unsigned ncont_interpolated_values() const
Number of continuously interpolated values:
unsigned Smoothed_vorticity_index
Index of smoothed vorticity – followed by derivatives; in 2D this has value 3.
unsigned npartial_derivative(const unsigned &n) const
Call the function written in VorticityRecoveryHelpers.
std::pair< unsigned, unsigned > recovered_dof_to_container_id(const unsigned &i) const
Helper function that, given the local dof number of the i-th vorticity or velocity derivative,...
void get_raw_vorticity_second_deriv(const Vector< double > &s, double &dvorticity_dxdy, const unsigned &index) const
Get raw derivative of smoothed derivative vorticity [0]: d^2/dx^2, [1]: d^2/dxdy, [2]: d^2/dy^2.
void vorticity_and_its_derivs(const Vector< double > &s, Vector< Vector< double > > &vort_and_derivs) const
Compute smoothed vorticity and its derivatives.
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 smoothed_vorticity_index() const
Index of smoothed vorticity – followed by derivatives.
void scalar_value_paraview(std::ofstream &file_out, const unsigned &i, const unsigned &nplot) const
Write values of the i-th scalar field at the plot points. Needs to be implemented for each new specif...
std::pair< unsigned, unsigned > vorticity_dof_to_container_id(const unsigned &i) const
Helper function that, given the local dof number of the i-th vorticity or velocity derivative,...
ExactVorticityFctPt Exact_vorticity_fct_pt
Pointer to function that specifies exact vorticity and derivs (for validation).
double vorticity_error_squared(const unsigned &i)
Compute the element's contribution to the (squared) L2 norm of the difference between exact and smoot...
unsigned Maximum_order_of_recoverable_velocity_derivatives
The current maximum order of velocity derivatives that can be recovered. Currently,...
unsigned get_maximum_order_of_recoverable_vorticity_derivative() const
The maximum order of vorticity derivative that can be recovered. This is set in the constructor and s...
unsigned Maximum_order_of_recoverable_vorticity_derivatives
The current maximum order of vorticity derivatives that can be recovered. Currently,...
unsigned nscalar_paraview() const
Number of scalars/fields output by this element. Re-implements broken virtual function in base class.
unsigned nvorticity_derivatives_to_recover() const
The number of terms calculated in the vorticity recovery. Also includes the zeroth derivative,...
void get_raw_velocity_deriv(const Vector< double > &s, Vector< double > &dveloc_dx) const
Get raw derivative of velocity.
void output(std::ostream &outfile, const unsigned &nplot)
Overloaded output function: Output velocity, pressure and the smoothed vorticity.
int Maximum_order_of_velocity_derivative
Maximum number of derivatives to retain in the velocity recovery. Note, the value 0 means we don't ca...
unsigned stored_dof_to_recoverable_dof(const unsigned &i) const
Given the STORED dof number, this function returns the global recovered number. For example,...
void output_smoothed_vorticity(std::ostream &outfile, const unsigned &nplot)
Output the velocity, smoothed vorticity and derivatives.
int get_maximum_order_of_vorticity_derivative() const
The maximum order of derivatives calculated in the vorticity recovery. Note, this value can only be s...
void get_raw_vorticity_third_deriv(const Vector< double > &s, Vector< double > &dvorticity_dxdxdy) const
Get raw derivative of smoothed derivative vorticity [0]: d^3/dx^3, [1]: d^3/dx^2dy,...
void get_raw_vorticity_deriv(const Vector< double > &s, double &dvorticity_dx, const unsigned &index) const
Get raw derivative of smoothed vorticity.
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...
unsigned Number_of_values_per_field
Number of values per field; how many of the following do we want: u,v,p,omega,d/dx,...
void get_raw_vorticity_deriv(const Vector< double > &s, Vector< double > &dvorticity_dx) const
Get raw derivative of smoothed vorticity.
unsigned get_maximum_order_of_recoverable_velocity_derivative() const
The maximum order of velocity derivative that can be recovered. This is set in the constructor and sh...
int Maximum_order_of_vorticity_derivative
Maximum number of derivatives to retain in the vorticity recovery. Note, the value -1 means we ONLY o...
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...
unsigned nvelocity_derivatives_to_recover() const
The number of derivatives calculated in the velocity recovery. This does NOT include the zeroth deriv...
void output_analytical_veloc_and_vorticity(std::ostream &outfile, const unsigned &nplot)
Output exact velocity, vorticity, derivatives and indicator based on functions specified by two funct...
Vector< Vector< double > > create_container_for_vorticity_and_derivatives() const
Helper function to create a container for the vorticity and its partial derivatives....
void pin_smoothed_vorticity()
Pin all smoothed vorticity quantities.
ExactVorticityFctPt & exact_vorticity_fct_pt()
Access function: Pointer to function that specifies exact vorticity and derivatives (for validation).
void(* ExactVorticityFctPt)(const Vector< double > &x, Vector< Vector< double > > &vort_and_derivs)
Typedef for pointer to function that specifies the exact vorticity and derivatives (for validation)
void get_raw_velocity_deriv(const Vector< double > &s, double &dveloc_dx, const unsigned &index) const
Get raw derivative of velocity.
void get_raw_vorticity_third_deriv(const Vector< double > &s, double &dvorticity_dxdxdy, const unsigned &index) const
Get raw derivative of smoothed derivative vorticity [0]: d^3/dx^3, [1]: d^3/dx^2dy,...
void get_raw_vorticity_second_deriv(const Vector< double > &s, Vector< double > &dvorticity_dxdy) const
Get raw derivative of smoothed derivative vorticity.
std::string scalar_name_paraview(const unsigned &i) const
Name of the i-th scalar field. Default implementation returns V1 for the first one,...
ExactVorticityFctPt exact_vorticity_fct_pt() const
Access function: Pointer to function that specifies exact vorticity and derivatives (for validation) ...
unsigned N_dim
Number of dimensions in the element.
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
unsigned & recovery_order()
Access function for order of recovery polynomials.
unsigned nrecovery_order() const
Get the recovery order.
void recover_vorticity(Mesh *mesh_pt, DocInfo &doc_info)
Recover vorticity from patches – output intermediate steps to directory specified by DocInfo object.
void recover_vorticity(Mesh *mesh_pt)
Recover vorticity from patches.
virtual ~VorticitySmoother()
Empty virtual destructor.
VorticitySmoother(const VorticitySmoother &)=delete
Broken copy constructor.
Integral * integral_rec(const bool &is_q_mesh)
Integation scheme associated with the recovery shape functions must be of sufficiently high order to ...
void shape_rec(const Vector< double > &x, Vector< double > &psi_r)
Recovery shape functions as functions of the global, Eulerian coordinate x of dimension dim....
VorticitySmoother(const unsigned &recovery_order)
Constructor: Set order of recovery shape functions.
void operator=(const VorticitySmoother &)=delete
Broken assignment operator.
unsigned Recovery_order
Order of recovery polynomials.
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)
Given the vector of elements that make up a patch, compute the vector of recovered vorticity coeffici...
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...
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
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:564
std::string to_string(T object, unsigned float_precision=8)
Conversion function that should work for anything with operator<< defined (at least all basic types).
double timer()
returns the time in seconds after some point in past
class oomph::VorticityRecoveryHelpers::RecoveryHelper Recovery_helper
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...