boussinesq_elements.h
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 //Header for a multi-physics elements that couple the Navier--Stokes
31 //and advection diffusion elements via a multi domain approach,
32 //giving Boussinesq convection
33 
34 #ifndef OOMPH_BOUSSINESQ_ELEMENTS_HEADER
35 #define OOMPH_BOUSSINESQ_ELEMENTS_HEADER
36 
37 //Oomph-lib headers, we require the generic, advection-diffusion
38 //and navier-stokes elements.
39 #include "generic.h"
40 #include "advection_diffusion.h"
41 #include "navier_stokes.h"
42 
43 
44 
45 namespace oomph
46 {
47 
48 /////////////////////////////////////////////////////////////////////////
49 /////////////////////////////////////////////////////////////////////////
50 /////////////////////////////////////////////////////////////////////////
51 
52 
53 //======================class definition==============================
54 ///A class that solves the Boussinesq approximation of the Navier--Stokes
55 ///and energy equations by coupling two pre-existing classes.
56 ///The QAdvectionDiffusionElement with bi-quadratic interpolation for the
57 ///scalar variable (temperature) and
58 ///QCrouzeixRaviartElement which solves the Navier--Stokes equations
59 ///using bi-quadratic interpolation for the velocities and a discontinuous
60 ///bi-linear interpolation for the pressure. Note that we are free to
61 ///choose the order in which we store the variables at the nodes. In this
62 ///case we choose to store the variables in the order fluid velocities
63 ///followed by temperature. We must, therefore, overload the function
64 ///AdvectionDiffusionEquations<DIM>::u_index_adv_diff() to indicate that
65 ///the temperature is stored at the DIM-th position not the 0-th. We do not
66 ///need to overload the corresponding function in the
67 ///NavierStokesEquations<DIM> class because the velocities are stored
68 ///first.
69 //=========================================================================
70 template<unsigned DIM>
72  public virtual QAdvectionDiffusionElement<DIM,3>,
73  public virtual QCrouzeixRaviartElement<DIM>
74 {
75 
76 private:
77 
78  /// Pointer to a private data member, the Rayleigh number
79  double* Ra_pt;
80 
81  /// The static default value of the Rayleigh number
83 
84 public:
85 
86  /// \short Constructor: call the underlying constructors and
87  /// initialise the pointer to the Rayleigh number to point
88  /// to the default value of 0.0.
91  {
93  }
94 
95  /// Unpin p_dof-th pressure dof
96  void unfix_pressure(const unsigned &p_dof)
97  {
98  this->internal_data_pt(this->P_nst_internal_index)->unpin(p_dof);
99  }
100 
101  ///\short The required number of values stored at the nodes is the sum of the
102  ///required values of the two single-physics elements. Note that this step is
103  ///generic for any multi-physics element of this type.
104  unsigned required_nvalue(const unsigned &n) const
107 
108  ///Access function for the Rayleigh number (const version)
109  const double &ra() const {return *Ra_pt;}
110 
111  ///Access function for the pointer to the Rayleigh number
112  double* &ra_pt() {return Ra_pt;}
113 
114  /// Final override for disable ALE
115  void disable_ALE()
116  {
117  //Disable ALE in both sets of equations
120  }
121 
122  /// Final override for enable ALE
123  void enable_ALE()
124  {
125  //Enable ALE in both sets of equations
128  }
129 
130  /// \short Number of scalars/fields output by this element. Broken
131  /// virtual. Needs to be implemented for each new specific element type.
132  /// Temporary dummy
133  unsigned nscalar_paraview() const
134  {
135  throw OomphLibError(
136  "This function hasn't been implemented for this element",
137  OOMPH_CURRENT_FUNCTION,
138  OOMPH_EXCEPTION_LOCATION);
139 
140  // Dummy unsigned
141  return 0;
142  }
143 
144  /// \short Write values of the i-th scalar field at the plot points. Broken
145  /// virtual. Needs to be implemented for each new specific element type.
146  /// Temporary dummy
147  void scalar_value_paraview(std::ofstream& file_out,
148  const unsigned& i,
149  const unsigned& nplot) const
150  {
151  throw OomphLibError(
152  "This function hasn't been implemented for this element",
153  OOMPH_CURRENT_FUNCTION,
154  OOMPH_EXCEPTION_LOCATION);
155  }
156 
157 
158  /// \short Name of the i-th scalar field. Default implementation
159  /// returns V1 for the first one, V2 for the second etc. Can (should!) be
160  /// overloaded with more meaningful names.
161  std::string scalar_name_paraview(const unsigned& i) const
162  {
163  return "V"+StringConversion::to_string(i);
164  }
165 
166 
167  /// Overload the standard output function with the broken default
168  void output(std::ostream &outfile) {FiniteElement::output(outfile);}
169 
170  /// \short Output function:
171  /// Output x, y, u, v, p, theta at Nplot^DIM plot points
172  // Start of output function
173  void output(std::ostream &outfile, const unsigned &nplot)
174  {
175  //vector of local coordinates
176  Vector<double> s(DIM);
177 
178  // Tecplot header info
179  outfile << this->tecplot_zone_string(nplot);
180 
181  // Loop over plot points
182  unsigned num_plot_points=this->nplot_points(nplot);
183  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
184  {
185  // Get local coordinates of plot point
186  this->get_s_plot(iplot,nplot,s);
187 
188  // Output the position of the plot point
189  for(unsigned i=0;i<DIM;i++)
190  {outfile << this->interpolated_x(s,i) << " ";}
191 
192  // Output the fluid velocities at the plot point
193  for(unsigned i=0;i<DIM;i++)
194  {outfile << this->interpolated_u_nst(s,i) << " ";}
195 
196  // Output the fluid pressure at the plot point
197  outfile << this->interpolated_p_nst(s) << " ";
198 
199  // Output the temperature (the advected variable) at the plot point
200  outfile << this->interpolated_u_adv_diff(s) << std::endl;
201  }
202  outfile << std::endl;
203 
204  // Write tecplot footer (e.g. FE connectivity lists)
205  this->write_tecplot_zone_footer(outfile,nplot);
206  } //End of output function
207 
208 
209  /// \short C-style output function: Broken default
210  void output(FILE* file_pt)
211  {FiniteElement::output(file_pt);}
212 
213  /// \short C-style output function: Broken default
214  void output(FILE* file_pt, const unsigned &n_plot)
215  {FiniteElement::output(file_pt,n_plot);}
216 
217  /// \short Output function for an exact solution: Broken default
218  void output_fct(std::ostream &outfile, const unsigned &Nplot,
220  exact_soln_pt)
221  {FiniteElement::output_fct(outfile,Nplot,exact_soln_pt);}
222 
223 
224  /// \short Output function for a time-dependent exact solution:
225  /// Broken default.
226  void output_fct(std::ostream &outfile, const unsigned &Nplot,
227  const double& time,
229  exact_soln_pt)
230  {
232  output_fct(outfile,Nplot,time,exact_soln_pt);
233  }
234 
235  ///\short Overload the index at which the temperature
236  ///variable is stored. We choose to store it after the fluid velocities.
237  inline unsigned u_index_adv_diff() const {return DIM;}
238 
239  /// \short Validate against exact solution at given time
240  /// Solution is provided via function pointer.
241  /// Plot at a given number of plot points and compute L2 error
242  /// and L2 norm of velocity solution over element
243  /// Call the broken default
244  void compute_error(std::ostream &outfile,
246  const double& time,
247  double& error, double& norm)
248  {FiniteElement::compute_error(outfile,exact_soln_pt,
249  time,error,norm);}
250 
251  /// \short Validate against exact solution.
252  /// Solution is provided via function pointer.
253  /// Plot at a given number of plot points and compute L2 error
254  /// and L2 norm of velocity solution over element
255  /// Call the broken default
256  void compute_error(std::ostream &outfile,
258  double& error, double& norm)
259  {FiniteElement::compute_error(outfile,exact_soln_pt,error,norm);}
260 
261  /// \short Overload the wind function in the advection-diffusion equations.
262  /// This provides the coupling from the Navier--Stokes equations to the
263  /// advection-diffusion equations because the wind is the fluid velocity.
264  void get_wind_adv_diff(const unsigned& ipt,
265  const Vector<double> &s, const Vector<double>& x,
266  Vector<double>& wind) const
267  {
268  //The wind function is simply the velocity at the points
269  this->interpolated_u_nst(s,wind);
270  }
271 
272 
273  /// \short Overload the body force in the Navier-Stokes equations
274  /// This provides the coupling from the advection-diffusion equations
275  /// to the Navier--Stokes equations, the body force is the
276  /// temperature multiplied by the Rayleigh number acting in the
277  /// direction opposite to gravity.
278  void get_body_force_nst(const double& time, const unsigned& ipt,
279  const Vector<double> &s, const Vector<double> &x,
280  Vector<double> &result)
281  {
282  //Get the temperature
283  const double interpolated_t = this->interpolated_u_adv_diff(s);
284 
285  // Get vector that indicates the direction of gravity from
286  // the Navier-Stokes equations
288 
289  // Temperature-dependent body force:
290  for (unsigned i=0;i<DIM;i++)
291  {
292  result[i] = -gravity[i]*interpolated_t*ra();
293  }
294  } // end of get_body_force
295 
296  /// \short Calculate the element's contribution to the residual vector.
297  /// Recall that fill_in_* functions MUST NOT initialise the entries
298  /// in the vector to zero. This allows us to call the
299  /// fill_in_* functions of the constituent single-physics elements
300  /// sequentially, without wiping out any previously computed entries.
302  {
303  //Fill in the residuals of the Navier-Stokes equations
305 
306  //Fill in the residuals of the advection-diffusion eqautions
308  residuals);
309  }
310 
311 
312 //-----------Finite-difference the entire jacobian-----------------------
313 //-----------------------------------------------------------------------
314 #ifdef USE_FD_JACOBIAN_FOR_BUOYANT_Q_ELEMENT
315 
316 
317  ///\short Compute the element's residual vector and the Jacobian matrix.
318  /// Jacobian is computed by finite-differencing.
320  DenseMatrix<double> &jacobian)
321  {
322  // This function computes the Jacobian by finite-differencing
324  }
325 
326  /// Add the element's contribution to its residuals vector,
327  /// jacobian matrix and mass matrix
329  Vector<double> &residuals, DenseMatrix<double> &jacobian,
330  DenseMatrix<double> &mass_matrix)
331  {
332  //Call the standard (Broken) function
333  //which will prevent these elements from being used
334  //in eigenproblems until replaced.
336  residuals,jacobian,mass_matrix);
337  }
338 
339 #else
340 //--------Finite-difference off-diagonal-blocks in jacobian--------------
341 //-----------------------------------------------------------------------
342 #ifdef USE_OFF_DIAGONAL_FD_JACOBIAN_FOR_BUOYANT_Q_ELEMENT
343 
344  ///\short Helper function to get the off-diagonal blocks of the Jacobian
345  ///matrix by finite differences
347  DenseMatrix<double> &jacobian)
348  {
349  //Local storage for the index in the nodes at which the
350  //Navier-Stokes velocities are stored (we know that this should be 0,1,2)
351  unsigned u_nodal_nst[DIM];
352  for(unsigned i=0;i<DIM;i++)
353  {u_nodal_nst[i] = this->u_index_nst(i);}
354 
355  //Local storage for the index at which the temperature is stored
356  unsigned u_nodal_adv_diff = this->u_index_adv_diff();
357 
358  //Find the total number of unknowns in the elements
359  unsigned n_dof = this->ndof();
360 
361  //Temporary storage for residuals
362  Vector<double> newres(n_dof);
363 
364  //Storage for local unknown and local equation
365  int local_unknown =0, local_eqn = 0;
366 
367  //Set the finite difference step
369 
370  //Find the number of nodes
371  unsigned n_node = this->nnode();
372 
373  //Calculate the contribution of the Navier--Stokes velocities to the
374  //advection-diffusion equations
375 
376  //Loop over the nodes
377  for(unsigned n=0;n<n_node;n++)
378  {
379  //There are DIM values of the velocities
380  for(unsigned i=0;i<DIM;i++)
381  {
382  //Get the local velocity equation number
383  local_unknown = this->nodal_local_eqn(n,u_nodal_nst[i]);
384 
385  //If it's not pinned
386  if(local_unknown >= 0)
387  {
388  //Get a pointer to the velocity value
389  double *value_pt = this->node_pt(n)->value_pt(u_nodal_nst[i]);
390 
391  //Save the old value
392  double old_var = *value_pt;
393 
394  //Increment the value
395  *value_pt += fd_step;
396 
397  //Get the altered advection-diffusion residuals.
398  //which must be done using fill_in_contribution because
399  //get_residuals will always return the full residuals
400  //because the appropriate fill_in function is overloaded above
401  for(unsigned m=0;m<n_dof;m++) {newres[m] = 0.0;}
404 
405  //AdvectionDiffusionEquations<DIM>::get_residuals(newres);
406 
407  //Now fill in the Advection-Diffusion sub-block
408  //of the jacobian
409  for(unsigned m=0;m<n_node;m++)
410  {
411  //Local equation for temperature
412  local_eqn = this->nodal_local_eqn(m,u_nodal_adv_diff);
413 
414  //If it's not a boundary condition
415  if(local_eqn >= 0)
416  {
417  double sum = (newres[local_eqn] - residuals[local_eqn])/fd_step;
418  jacobian(local_eqn,local_unknown) = sum;
419  }
420  }
421 
422  //Reset the nodal data
423  *value_pt = old_var;
424  }
425  }
426 
427 
428  //Calculate the contribution of the temperature to the Navier--Stokes
429  //equations
430  {
431  //Get the local equation number
432  local_unknown = this->nodal_local_eqn(n,u_nodal_adv_diff);
433 
434  //If it's not pinned
435  if(local_unknown >= 0)
436  {
437  //Get a pointer to the concentration value
438  double *value_pt = this->node_pt(n)->value_pt(u_nodal_adv_diff);
439 
440  //Save the old value
441  double old_var = *value_pt;
442 
443  //Increment the value (Really need access)
444  *value_pt += fd_step;
445 
446  //Get the altered Navier--Stokes residuals
447  //which must be done using fill_in_contribution because
448  //get_residuals will always return the full residuals
449  //because the appropriate fill_in function is overloaded above
450  for(unsigned m=0;m<n_dof;m++) {newres[m] = 0.0;}
453 
454  //NavierStokesEquations<DIM>::get_residuals(newres);
455 
456  //Now fill in the Navier-Stokes sub-block
457  for(unsigned m=0;m<n_node;m++)
458  {
459  //Loop over the fluid velocities
460  for(unsigned j=0;j<DIM;j++)
461  {
462  //Local fluid equation
463  local_eqn = this->nodal_local_eqn(m,u_nodal_nst[j]);
464  if(local_eqn >= 0)
465  {
466  double sum = (newres[local_eqn] - residuals[local_eqn])/fd_step;
467  jacobian(local_eqn,local_unknown) = sum;
468  }
469  }
470  }
471 
472  //Reset the nodal data
473  *value_pt = old_var;
474  }
475  }
476 
477  } //End of loop over nodes
478  }
479 
480 
481  ///\short Compute the element's residual Vector and the Jacobian matrix.
482  /// Use finite-differencing only for the off-diagonal blocks.
484  DenseMatrix<double> &jacobian)
485  {
486 
487  //Calculate the Navier-Stokes contributions (diagonal block and residuals)
489  fill_in_contribution_to_jacobian(residuals,jacobian);
490 
491  //Calculate the advection-diffusion contributions
492  //(diagonal block and residuals)
494  fill_in_contribution_to_jacobian(residuals,jacobian);
495 
496  //We now fill in the off-diagonal (interaction) blocks by finite
497  //differences.
498  fill_in_off_diagonal_jacobian_blocks_by_fd(residuals,jacobian);
499  } //End of jacobian calculation
500 
501 
502  /// Add the element's contribution to its residuals vector,
503  /// jacobian matrix and mass matrix
505  Vector<double> &residuals, DenseMatrix<double> &jacobian,
506  DenseMatrix<double> &mass_matrix)
507  {
508  //Get the analytic diagonal terms
511  residuals,jacobian,mass_matrix);
512 
515  residuals,jacobian,mass_matrix);
516 
517  //Now fill in the off-diagonal blocks
518  fill_in_off_diagonal_jacobian_blocks_by_fd(residuals,jacobian);
519  }
520 
521 
522  //--------------------Analytic jacobian---------------------------------
523 //-----------------------------------------------------------------------
524 #else
525 
526  ///\short Helper function to get the off-diagonal blocks of the Jacobian
527  ///matrix analytically
529  Vector<double> &residuals, DenseMatrix<double> &jacobian)
530  {
531  //We now fill in the off-diagonal (interaction) blocks analytically
532  //This requires knowledge of exactly how the residuals are assembled
533  //within the parent elements and involves yet another loop over
534  //the integration points!
535 
536  //Local storage for the index in the nodes at which the
537  //Navier-Stokes velocities are stored (we know that this should be 0,1,2)
538  unsigned u_nodal_nst[DIM];
539  for(unsigned i=0;i<DIM;i++) {u_nodal_nst[i] = this->u_index_nst(i);}
540 
541  //Local storage for the index at which the temperature is stored
542  const unsigned u_nodal_adv_diff = this->u_index_adv_diff();
543 
544  //Find out how many nodes there are
545  const unsigned n_node = this->nnode();
546 
547  //Set up memory for the shape and test functions and their derivatives
548  Shape psif(n_node), testf(n_node);
549  DShape dpsifdx(n_node,DIM), dtestfdx(n_node,DIM);
550 
551  //Number of integration points
552  const unsigned n_intpt = this->integral_pt()->nweight();
553 
554  //Get Physical Variables from Element
555  double Ra = this->ra();
556  double Pe = this->pe();
557  Vector<double> gravity = this->g();
558 
559  //Integers to store the local equations and unknowns
560  int local_eqn=0, local_unknown=0;
561 
562  //Loop over the integration points
563  for(unsigned ipt=0;ipt<n_intpt;ipt++)
564  {
565  //Get the integral weight
566  double w = this->integral_pt()->weight(ipt);
567 
568  //Call the derivatives of the shape and test functions
569  double J =
570  this->dshape_and_dtest_eulerian_at_knot_nst(ipt,psif,dpsifdx,
571  testf,dtestfdx);
572 
573  //Premultiply the weights and the Jacobian
574  double W = w*J;
575 
576  //Calculate local values of temperature derivatives
577  //Allocate
578  Vector<double> interpolated_du_adv_diff_dx(DIM,0.0);
579 
580  // Loop over nodes
581  for(unsigned l=0;l<n_node;l++)
582  {
583  //Get the nodal value
584  double u_value = this->raw_nodal_value(l,u_nodal_adv_diff);
585  //Loop over the derivative directions
586  for(unsigned j=0;j<DIM;j++)
587  {
588  interpolated_du_adv_diff_dx[j] += u_value*dpsifdx(l,j);
589  }
590  }
591 
592  //Assemble the jacobian terms
593 
594  //Loop over the test functions
595  for(unsigned l=0;l<n_node;l++)
596  {
597  //Assemble the contributions of the temperature to
598  //the Navier--Stokes equations (which arise through the buoyancy
599  //body-force term)
600 
601  //Loop over the velocity components in the Navier--Stokes equtions
602  for(unsigned i=0;i<DIM;i++)
603  {
604  //If it's not a boundary condition
605  local_eqn = this->nodal_local_eqn(l,u_nodal_nst[i]);
606  if(local_eqn >= 0)
607  {
608  //Loop over the velocity shape functions again
609  for(unsigned l2=0;l2<n_node;l2++)
610  {
611  //We have only the temperature degree of freedom to consider
612  //If it's non-zero add in the contribution
613  local_unknown = this->nodal_local_eqn(l2,u_nodal_adv_diff);
614  if(local_unknown >= 0)
615  {
616  //Add contribution to jacobian matrix
617  jacobian(local_eqn,local_unknown)
618  += -gravity[i]*psif(l2)*Ra*testf(l)*W;
619  }
620  }
621  }
622  }
623 
624  //Assemble the contributions of the fluid velocities to the
625  //advection-diffusion equation for the temperature
626  {
627  local_eqn = this->nodal_local_eqn(l,u_nodal_adv_diff);
628  //If it's not pinned
629  if(local_eqn >= 0)
630  {
631  //Loop over the shape functions again
632  for(unsigned l2=0;l2<n_node;l2++)
633  {
634  //Loop over the velocity degrees of freedom
635  for(unsigned i2=0;i2<DIM;i2++)
636  {
637  //Get the local unknown
638  local_unknown = this->nodal_local_eqn(l2,u_nodal_nst[i2]);
639  //If it's not pinned
640  if(local_unknown >= 0)
641  {
642  //Add the contribution to the jacobian matrix
643  jacobian(local_eqn,local_unknown)
644  -= Pe*psif(l2)*interpolated_du_adv_diff_dx[i2]*testf(l)*W;
645  }
646  }
647  }
648  }
649  }
650 
651  }
652  }
653  }
654 
655 
656  ///\short Compute the element's residual Vector and the Jacobian matrix.
657  /// Use analytic expressions for the off-diagonal blocks
659  DenseMatrix<double> &jacobian)
660  {
661 
662  //Calculate the Navier-Stokes contributions (diagonal block and residuals)
664  fill_in_contribution_to_jacobian(residuals,jacobian);
665 
666  //Calculate the advection-diffusion contributions
667  //(diagonal block and residuals)
669  fill_in_contribution_to_jacobian(residuals,jacobian);
670 
671  //Fill in the off diagonal blocks analytically
673  }
674 
675 
676  /// Add the element's contribution to its residuals vector,
677  /// jacobian matrix and mass matrix
679  Vector<double> &residuals, DenseMatrix<double> &jacobian,
680  DenseMatrix<double> &mass_matrix)
681  {
682  //Get the analytic diagonal terms for the jacobian and mass matrix
685  residuals,jacobian,mass_matrix);
686 
689  residuals,jacobian,mass_matrix);
690 
691  //Now fill in the off-diagonal blocks in the jacobian matrix
693  }
694 
695 #endif
696 #endif
697 
698 };
699 
700 //=========================================================================
701 /// Set the default physical value to be zero
702 //=========================================================================
703 template<>
705 template<>
707 
708 
709 //=======================================================================
710 /// Face geometry of the 2D Buoyant Crouzeix_Raviart elements
711 //=======================================================================
712 template<unsigned int DIM>
714 public virtual QElement<DIM-1,3>
715 {
716  public:
717  FaceGeometry() : QElement<DIM-1,3>() {}
718 };
719 
720 //=======================================================================
721 /// Face geometry of the Face geometry of 2D Buoyant Crouzeix_Raviart elements
722 //=======================================================================
723 template<>
725 public virtual PointElement
726 {
727  public:
729 };
730 
731 
732 /////////////////////////////////////////////////////////////////////////
733 /////////////////////////////////////////////////////////////////////////
734 /////////////////////////////////////////////////////////////////////////
735 
736 
737 
738 //============start_element_class============================================
739 ///A RefineableElement class that solves the
740 ///Boussinesq approximation of the Navier--Stokes
741 ///and energy equations by coupling two pre-existing classes.
742 ///The RefineableQAdvectionDiffusionElement
743 ///with bi-quadratic interpolation for the
744 ///scalar variable (temperature) and
745 ///RefineableQCrouzeixRaviartElement which solves the Navier--Stokes equations
746 ///using bi-quadratic interpolation for the velocities and a discontinuous
747 ///bi-linear interpolation for the pressure. Note that we are free to
748 ///choose the order in which we store the variables at the nodes. In this
749 ///case we choose to store the variables in the order fluid velocities
750 ///followed by temperature. We must, therefore, overload the function
751 ///AdvectionDiffusionEquations<DIM>::u_index_adv_diff() to indicate that
752 ///the temperature is stored at the DIM-th position not the 0-th. We do not
753 ///need to overload the corresponding function in the
754 ///NavierStokesEquations<DIM> class because the velocities are stored
755 ///first. Finally, we choose to use the flux-recovery calculation from the
756 ///fluid velocities to provide the error used in the mesh adaptation.
757 //==========================================================================
758 template<unsigned DIM>
760 public virtual RefineableQAdvectionDiffusionElement<DIM,3>,
761 public virtual RefineableQCrouzeixRaviartElement<DIM>
762 {
763 
764 private:
765 
766  ///Pointer to a new physical variable, the Rayleigh number
767  double* Ra_pt;
768 
769  /// The static default value of the Rayleigh number
771 
772 public:
773  ///\short Constructor: call the underlying constructors and
774  ///initialise the pointer to the Rayleigh number to address the default
775  ///value of 0.0
779  {
781  }
782 
783  ///\short The required number of values stored at the nodes is
784  ///the sum of the required values of the two single-physics elements. This
785  ///step is generic for any composed element of this type.
786  inline unsigned required_nvalue(const unsigned &n) const
789 
790  ///Access function for the Rayleigh number (const version)
791  const double &ra() const {return *Ra_pt;}
792 
793  ///Access function for the pointer to the Rayleigh number
794  double* &ra_pt() {return Ra_pt;}
795 
796 
797  /// Final override for disable ALE
798  void disable_ALE()
799  {
800  //Disable ALE in both sets of equations
803  }
804 
805  /// Final override for enable ALE
806  void enable_ALE()
807  {
808  //Enable ALE in both sets of equations
811  }
812 
813 
814  /// \short Number of scalars/fields output by this element. Broken
815  /// virtual. Needs to be implemented for each new specific element type.
816  /// Temporary dummy
817  unsigned nscalar_paraview() const
818  {
819  throw OomphLibError(
820  "This function hasn't been implemented for this element",
821  OOMPH_CURRENT_FUNCTION,
822  OOMPH_EXCEPTION_LOCATION);
823 
824  // Dummy unsigned
825  return 0;
826  }
827 
828  /// \short Write values of the i-th scalar field at the plot points. Broken
829  /// virtual. Needs to be implemented for each new specific element type.
830  /// Temporary dummy
831  void scalar_value_paraview(std::ofstream& file_out,
832  const unsigned& i,
833  const unsigned& nplot) const
834  {
835  throw OomphLibError(
836  "This function hasn't been implemented for this element",
837  OOMPH_CURRENT_FUNCTION,
838  OOMPH_EXCEPTION_LOCATION);
839  }
840 
841 
842  /// \short Name of the i-th scalar field. Default implementation
843  /// returns V1 for the first one, V2 for the second etc. Can (should!) be
844  /// overloaded with more meaningful names.
845  std::string scalar_name_paraview(const unsigned& i) const
846  {
847  return "V"+StringConversion::to_string(i);
848  }
849 
850  /// Overload the standard output function with the broken default
851  void output(std::ostream &outfile)
852  {FiniteElement::output(outfile);}
853 
854  /// \short Output function:
855  /// x,y,u or x,y,z,u at Nplot^DIM plot points
856  void output(std::ostream &outfile, const unsigned &nplot)
857  {
858  //vector of local coordinates
859  Vector<double> s(DIM);
860 
861  // Tecplot header info
862  outfile << this->tecplot_zone_string(nplot);
863 
864  // Loop over plot points
865  unsigned num_plot_points=this->nplot_points(nplot);
866  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
867  {
868  // Get local coordinates of plot point
869  this->get_s_plot(iplot,nplot,s);
870 
871  // Output the position of the plot point
872  for(unsigned i=0;i<DIM;i++)
873  {outfile << this->interpolated_x(s,i) << " ";}
874 
875  // Output the fluid velocities at the plot point
876  for(unsigned i=0;i<DIM;i++)
877  {outfile << this->interpolated_u_nst(s,i) << " ";}
878 
879  // Output the fluid pressure at the plot point
880  outfile << this->interpolated_p_nst(s) << " ";
881 
882  // Output the temperature at the plot point
883  outfile << this->interpolated_u_adv_diff(s) << " " << std::endl;
884  }
885  outfile << std::endl;
886 
887  // Write tecplot footer (e.g. FE connectivity lists)
888  this->write_tecplot_zone_footer(outfile,nplot);
889 
890  }
891 
892  /// \short C-style output function: Broken default
893  void output(FILE* file_pt)
894  {FiniteElement::output(file_pt);}
895 
896  /// \short C-style output function: Broken default
897  void output(FILE* file_pt, const unsigned &n_plot)
898  {FiniteElement::output(file_pt,n_plot);}
899 
900  /// \short Output function for an exact solution: Broken default
901  void output_fct(std::ostream &outfile, const unsigned &Nplot,
903  exact_soln_pt)
904  {FiniteElement::output_fct(outfile,Nplot,exact_soln_pt);}
905 
906 
907  /// \short Output function for a time-dependent exact solution.
908  /// Broken default
909  void output_fct(std::ostream &outfile, const unsigned &Nplot,
910  const double& time,
912  exact_soln_pt)
913  {
915  output_fct(outfile,Nplot,time,exact_soln_pt);
916  }
917 
918  ///\short Overload the index at which the temperature
919  ///variable is stored. We choose to store is after the fluid velocities.
920  unsigned u_index_adv_diff() const
921  {return DIM;}
922 
923  /// \short Number of vertex nodes in the element is obtained from the
924  /// geometric element.
925  unsigned nvertex_node() const
927 
928  /// \short Pointer to the j-th vertex node in the element,
929  /// Call the geometric element's function.
930  Node* vertex_node_pt(const unsigned& j) const
932 
933  /// \short The total number of continously interpolated values is
934  /// DIM+1 (DIM fluid velocities and one temperature).
935  unsigned ncont_interpolated_values() const
936  {return DIM+1;}
937 
938 
939  /// \short Get the continuously interpolated values at the local coordinate s.
940  /// We choose to put the fluid velocities first, followed by the
941  /// temperature.
943  {
944  //Storage for the fluid velocities
945  Vector<double> nst_values;
946 
947  //Get the fluid velocities from the fluid element
949  get_interpolated_values(s,nst_values);
950 
951  //Storage for the temperature
952  Vector<double> advection_values;
953 
954  //Get the temperature from the advection-diffusion element
956  get_interpolated_values(s,advection_values);
957 
958  //Add the fluid velocities to the values vector
959  for(unsigned i=0;i<DIM;i++) {values.push_back(nst_values[i]);}
960 
961  //Add the concentration to the end
962  values.push_back(advection_values[0]);
963  }
964 
965 
966 
967  /// \short Get all continuously interpolated values at the local
968  /// coordinate s at time level t (t=0: present; t>0: previous).
969  /// We choose to put the fluid velocities first, followed by the
970  /// temperature
971  void get_interpolated_values(const unsigned& t, const Vector<double>&s,
972  Vector<double>& values)
973  {
974  //Storage for the fluid velocities
975  Vector<double> nst_values;
976 
977  //Get the fluid velocities from the fluid element
979  get_interpolated_values(t,s,nst_values);
980 
981  //Storage for the temperature
982  Vector<double> advection_values;
983 
984  //Get the temperature from the advection-diffusion element
986  get_interpolated_values(s,advection_values);
987 
988  //Add the fluid velocities to the values vector
989  for(unsigned i=0;i<DIM;i++) {values.push_back(nst_values[i]);}
990 
991  //Add the concentration to the end
992  values.push_back(advection_values[0]);
993 
994  } // end of get_interpolated_values
995 
996 
997 
998  /// \short The additional hanging node information must be set up
999  /// for both single-physics elements.
1001  {
1004  }
1005 
1006 
1007 
1008  /// \short Call the rebuild_from_sons functions for each of the
1009  /// constituent multi-physics elements.
1010  void rebuild_from_sons(Mesh* &mesh_pt)
1011  {
1014  }
1015 
1016 
1017 
1018  /// \short Call the underlying single-physics element's further_build()
1019  /// functions and make sure that the pointer to the Rayleigh number
1020  /// is passed to the sons
1022  {
1025 
1026  //Cast the pointer to the father element to the specific
1027  //element type
1028  RefineableBuoyantQCrouzeixRaviartElement<DIM>* cast_father_element_pt
1030  this->father_element_pt());
1031 
1032  //Set the pointer to the Rayleigh number to be the same as that in
1033  //the father
1034  this->Ra_pt = cast_father_element_pt->ra_pt();
1035  } //end of further build
1036 
1037 
1038 
1039  /// The recovery order is that of the NavierStokes elements.
1040  unsigned nrecovery_order()
1042 
1043  /// \short The number of Z2 flux terms is the same as that in
1044  /// the fluid element plus that in the advection-diffusion element
1046  {
1049  }
1050 
1051 
1052  /// \short Get the Z2 flux by concatenating the fluxes from the fluid and
1053  /// the advection diffusion elements.
1055  {
1056  //Find the number of fluid fluxes
1057  unsigned n_fluid_flux =
1059 
1060  //Fill in the first flux entries as the velocity entries
1062 
1063  //Find the number of temperature fluxes
1064  unsigned n_temp_flux =
1066  Vector<double> temp_flux(n_temp_flux);
1067 
1068  //Get the temperature flux
1070 
1071  //Add the temperature flux to the end of the flux vector
1072  for(unsigned i=0;i<n_temp_flux;i++)
1073  {
1074  flux[n_fluid_flux+i] = temp_flux[i];
1075  }
1076 
1077  } //end of get_Z2_flux
1078 
1079  /// \short The number of compound fluxes is two (one for the fluid and
1080  /// one for the temperature)
1081  unsigned ncompound_fluxes() {return 2;}
1082 
1083  /// \short Fill in which flux components are associated with the fluid
1084  /// measure and which are associated with the temperature measure
1086  {
1087  //Find the number of fluid fluxes
1088  unsigned n_fluid_flux =
1090  //Find the number of temperature fluxes
1091  unsigned n_temp_flux =
1093 
1094  //The fluid fluxes are first
1095  //The values of the flux_index vector are zero on entry, so we
1096  //could omit this line
1097  for(unsigned i=0;i<n_fluid_flux;i++) {flux_index[i] = 0;}
1098 
1099  //Set the temperature fluxes (the last set of fluxes
1100  for(unsigned i=0;i<n_temp_flux;i++) {flux_index[n_fluid_flux + i] = 1;}
1101 
1102  } //end of get_Z2_compound_flux_indices
1103 
1104 
1105  /// \short Validate against exact solution at given time
1106  /// Solution is provided via function pointer.
1107  /// Plot at a given number of plot points and compute L2 error
1108  /// and L2 norm of velocity solution over element
1109  /// Overload to broken default
1110  void compute_error(std::ostream &outfile,
1112  const double& time,
1113  double& error, double& norm)
1114  {FiniteElement::compute_error(outfile,exact_soln_pt,
1115  time,error,norm);}
1116 
1117  /// \short Validate against exact solution.
1118  /// Solution is provided via function pointer.
1119  /// Plot at a given number of plot points and compute L2 error
1120  /// and L2 norm of velocity solution over element
1121  /// Overload to broken default.
1122  void compute_error(std::ostream &outfile,
1124  double& error, double& norm)
1126  compute_error(outfile,exact_soln_pt,error,norm);}
1127 
1128  /// \short Overload the wind function in the advection-diffusion equations.
1129  /// This provides the coupling from the Navier--Stokes equations to the
1130  /// advection-diffusion equations because the wind is the fluid velocity.
1131  void get_wind_adv_diff(const unsigned& ipt,
1132  const Vector<double> &s, const Vector<double>& x,
1133  Vector<double>& wind) const
1134  {
1135  //The wind function is simply the velocity at the points
1136  this->interpolated_u_nst(s,wind);
1137  }
1138 
1139 
1140  /// \short Overload the body force in the navier-stokes equations
1141  /// This provides the coupling from the advection-diffusion equations
1142  /// to the Navier--Stokes equations, the body force is the
1143  /// temperature multiplied by the Rayleigh number acting in the
1144  /// direction opposite to gravity.
1145  void get_body_force_nst(const double& time, const unsigned& ipt,
1146  const Vector<double> &s, const Vector<double> &x,
1147  Vector<double> &result)
1148  {
1149  // Get the temperature
1150  const double interpolated_t = this->interpolated_u_adv_diff(s);
1151 
1152  // Get vector that indicates the direction of gravity from
1153  // the Navier-Stokes equations
1155 
1156  // Temperature-dependent body force:
1157  for (unsigned i=0;i<DIM;i++)
1158  {
1159  result[i] = -gravity[i]*interpolated_t*ra();
1160  }
1161  }
1162 
1163  /// Fill in the constituent elements' contribution to the residual vector.
1165  {
1166  //Call the residuals of the Navier-Stokes equations
1168  residuals);
1169 
1170  //Call the residuals of the advection-diffusion equations
1173  }
1174 
1175 
1176  ///\short Compute the element's residual Vector and the jacobian matrix
1177  ///using full finite differences, the default implementation
1179  DenseMatrix<double> &jacobian)
1180  {
1181 #ifdef USE_FD_JACOBIAN_FOR_REFINEABLE_BUOYANT_Q_ELEMENT
1183 #else
1184  //Calculate the Navier-Stokes contributions (diagonal block and residuals)
1186  fill_in_contribution_to_jacobian(residuals,jacobian);
1187 
1188  //Calculate the advection-diffusion contributions
1189  //(diagonal block and residuals)
1191  fill_in_contribution_to_jacobian(residuals,jacobian);
1192 
1193  //We now fill in the off-diagonal (interaction) blocks analytically
1194  this->fill_in_off_diagonal_jacobian_blocks_analytic(residuals,jacobian);
1195 #endif
1196  } //End of jacobian calculation
1197 
1198  /// Add the element's contribution to its residuals vector,
1199  /// jacobian matrix and mass matrix
1201  Vector<double> &residuals, DenseMatrix<double> &jacobian,
1202  DenseMatrix<double> &mass_matrix)
1203  {
1204  //Call the (broken) version in the base class
1206  residuals,jacobian,mass_matrix);
1207  }
1208 
1209  /// \short Compute the contribution of the off-diagonal blocks
1210  /// analytically.
1212  Vector<double> &residuals, DenseMatrix<double> &jacobian)
1213  {
1214  //Perform another loop over the integration loops using the information
1215  //from the original elements' residual assembly loops to determine
1216  //the conributions to the jacobian
1217 
1218  // Local storage for pointers to hang_info objects
1219  HangInfo *hang_info_pt=0, *hang_info2_pt=0;
1220 
1221  //Local storage for the index in the nodes at which the
1222  //Navier-Stokes velocities are stored (we know that this should be 0,1,2)
1223  unsigned u_nodal_nst[DIM];
1224  for(unsigned i=0;i<DIM;i++) {u_nodal_nst[i] = this->u_index_nst(i);}
1225 
1226  //Local storage for the index at which the temperature is stored
1227  const unsigned u_nodal_adv_diff = this->u_index_adv_diff();
1228 
1229  //Find out how many nodes there are
1230  const unsigned n_node = this->nnode();
1231 
1232  //Set up memory for the shape and test functions and their derivatives
1233  Shape psif(n_node), testf(n_node);
1234  DShape dpsifdx(n_node,DIM), dtestfdx(n_node,DIM);
1235 
1236  //Number of integration points
1237  const unsigned n_intpt = this->integral_pt()->nweight();
1238 
1239  //Get Physical Variables from Element
1240  double Ra = this->ra();
1241  double Pe = this->pe();
1242  Vector<double> gravity = this->g();
1243 
1244  //Integers to store the local equations and unknowns
1245  int local_eqn=0, local_unknown=0;
1246 
1247  //Loop over the integration points
1248  for(unsigned ipt=0;ipt<n_intpt;ipt++)
1249  {
1250  //Get the integral weight
1251  double w = this->integral_pt()->weight(ipt);
1252 
1253  //Call the derivatives of the shape and test functions
1254  double J =
1255  this->dshape_and_dtest_eulerian_at_knot_nst(ipt,psif,dpsifdx,
1256  testf,dtestfdx);
1257 
1258  //Premultiply the weights and the Jacobian
1259  double W = w*J;
1260 
1261  //Calculate local values of temperature derivatives
1262  //Allocate
1263  Vector<double> interpolated_du_adv_diff_dx(DIM,0.0);
1264 
1265  // Loop over nodes
1266  for(unsigned l=0;l<n_node;l++)
1267  {
1268  //Get the nodal value
1269  double u_value = this->nodal_value(l,u_nodal_adv_diff);
1270  //Loop over the derivative directions
1271  for(unsigned j=0;j<DIM;j++)
1272  {
1273  interpolated_du_adv_diff_dx[j] += u_value*dpsifdx(l,j);
1274  }
1275  }
1276 
1277  //Assemble the Jacobian terms
1278  //---------------------------
1279 
1280  //Loop over the test functions/eqns
1281  for(unsigned l=0;l<n_node;l++)
1282  {
1283  //Local variables to store the number of master nodes and
1284  //the weight associated with the shape function if the node is hanging
1285  unsigned n_master=1;
1286  double hang_weight=1.0;
1287 
1288  //Local bool (is the node hanging)
1289  bool is_node_hanging = this->node_pt(l)->is_hanging();
1290 
1291  //If the node is hanging, get the number of master nodes
1292  if(is_node_hanging)
1293  {
1294  hang_info_pt = this->node_pt(l)->hanging_pt();
1295  n_master = hang_info_pt->nmaster();
1296  }
1297  //Otherwise there is just one master node, the node itself
1298  else
1299  {
1300  n_master = 1;
1301  }
1302 
1303  //Loop over the master nodes
1304  for(unsigned m=0;m<n_master;m++)
1305  {
1306  //If the node is hanging get weight from master node
1307  if(is_node_hanging)
1308  {
1309  //Get the hang weight from the master node
1310  hang_weight = hang_info_pt->master_weight(m);
1311  }
1312  else
1313  {
1314  // Node contributes with full weight
1315  hang_weight = 1.0;
1316  }
1317 
1318 
1319  //Assemble derivatives of Navier Stokes momentum w.r.t. temperature
1320  //-----------------------------------------------------------------
1321 
1322  // Loop over velocity components for equations
1323  for(unsigned i=0;i<DIM;i++)
1324  {
1325 
1326  //Get the equation number
1327  if(is_node_hanging)
1328  {
1329  //Get the equation number from the master node
1330  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
1331  u_nodal_nst[i]);
1332  }
1333  else
1334  {
1335  // Local equation number
1336  local_eqn = this->nodal_local_eqn(l,u_nodal_nst[i]);
1337  }
1338 
1339  if(local_eqn >= 0)
1340  {
1341  //Local variables to store the number of master nodes
1342  //and the weights associated with each hanging node
1343  unsigned n_master2=1;
1344  double hang_weight2=1.0;
1345 
1346  //Loop over the nodes for the unknowns
1347  for(unsigned l2=0;l2<n_node;l2++)
1348  {
1349  //Local bool (is the node hanging)
1350  bool is_node2_hanging = this->node_pt(l2)->is_hanging();
1351 
1352  //If the node is hanging, get the number of master nodes
1353  if(is_node2_hanging)
1354  {
1355  hang_info2_pt = this->node_pt(l2)->hanging_pt();
1356  n_master2 = hang_info2_pt->nmaster();
1357  }
1358  //Otherwise there is one master node, the node itself
1359  else
1360  {
1361  n_master2 = 1;
1362  }
1363 
1364  //Loop over the master nodes
1365  for(unsigned m2=0;m2<n_master2;m2++)
1366  {
1367  if(is_node2_hanging)
1368  {
1369  //Read out the local unknown from the master node
1370  local_unknown =
1371  this->local_hang_eqn(hang_info2_pt->master_node_pt(m2),
1372  u_nodal_adv_diff);
1373  //Read out the hanging weight from the master node
1374  hang_weight2 = hang_info2_pt->master_weight(m2);
1375  }
1376  else
1377  {
1378  //The local unknown number comes from the node
1379  local_unknown = this->nodal_local_eqn(l2,u_nodal_adv_diff);
1380  //The hang weight is one
1381  hang_weight2 = 1.0;
1382  }
1383 
1384  if(local_unknown >= 0)
1385  {
1386  //Add contribution to jacobian matrix
1387  jacobian(local_eqn,local_unknown)
1388  += -gravity[i]*psif(l2)*Ra*testf(l)*
1389  W*hang_weight*hang_weight2;
1390  }
1391  }
1392  }
1393  }
1394  }
1395 
1396 
1397  //Assemble derivative of adv diff eqn w.r.t. fluid veloc
1398  //------------------------------------------------------
1399  {
1400  //Get the equation number
1401  if(is_node_hanging)
1402  {
1403  //Get the equation number from the master node
1404  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
1405  u_nodal_adv_diff);
1406  }
1407  else
1408  {
1409  // Local equation number
1410  local_eqn = this->nodal_local_eqn(l,u_nodal_adv_diff);
1411  }
1412 
1413  //If it's not pinned
1414  if(local_eqn >= 0)
1415  {
1416  //Local variables to store the number of master nodes
1417  //and the weights associated with each hanging node
1418  unsigned n_master2=1;
1419  double hang_weight2=1.0;
1420 
1421  //Loop over the nodes for the unknowns
1422  for(unsigned l2=0;l2<n_node;l2++)
1423  {
1424  //Local bool (is the node hanging)
1425  bool is_node2_hanging = this->node_pt(l2)->is_hanging();
1426 
1427  //If the node is hanging, get the number of master nodes
1428  if(is_node2_hanging)
1429  {
1430  hang_info2_pt = this->node_pt(l2)->hanging_pt();
1431  n_master2 = hang_info2_pt->nmaster();
1432  }
1433  //Otherwise there is one master node, the node itself
1434  else
1435  {
1436  n_master2 = 1;
1437  }
1438 
1439  //Loop over the master nodes
1440  for(unsigned m2=0;m2<n_master2;m2++)
1441  {
1442  //If the node is hanging
1443  if(is_node2_hanging)
1444  {
1445  //Read out the hanging weight from the master node
1446  hang_weight2 = hang_info2_pt->master_weight(m2);
1447  }
1448  //If the node is not hanging
1449  else
1450  {
1451  //The hang weight is one
1452  hang_weight2 = 1.0;
1453  }
1454 
1455  //Loop over the velocity degrees of freedom
1456  for(unsigned i2=0;i2<DIM;i2++)
1457  {
1458  //If the node is hanging
1459  if(is_node2_hanging)
1460  {
1461  //Read out the local unknown from the master node
1462  local_unknown =
1463  this->local_hang_eqn(hang_info2_pt->master_node_pt(m2),
1464  u_nodal_nst[i2]);
1465  }
1466  else
1467  {
1468  //The local unknown number comes from the node
1469  local_unknown=this->nodal_local_eqn(l2, u_nodal_nst[i2]);
1470  }
1471 
1472  //If it's not pinned
1473  if(local_unknown >= 0)
1474  {
1475  //Add the contribution to the jacobian matrix
1476  jacobian(local_eqn,local_unknown)
1477  -= Pe*psif(l2)*interpolated_du_adv_diff_dx[i2]*testf(l)
1478  *W*hang_weight*hang_weight2;
1479  }
1480  }
1481  }
1482  }
1483  }
1484  }
1485 
1486  }
1487  }
1488  }
1489  } //End of function
1490 
1491 
1492 };
1493 
1494 
1495 //===================================================================
1496 ///Set the default value of the Rayleigh number to be zero
1497 //===================================================================
1498 template<>
1499 double RefineableBuoyantQCrouzeixRaviartElement<2>::
1501 
1502 template<>
1505 
1506 
1507 
1508 }
1509 
1510 #endif
void output_fct(std::ostream &outfile, const unsigned &Nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output function for an exact solution: Broken default.
unsigned num_Z2_flux_terms()
Number of 'flux' terms for Z2 error estimation.
int local_hang_eqn(Node *const &node_pt, const unsigned &i)
Access function that returns the local equation number for the hanging node variables (values stored ...
void unpin(const unsigned &i)
Unpin the i-th stored variable.
Definition: nodes.h:386
void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Get 'flux' for Z2 error recovery: Upper triangular entries in strain rate tensor. ...
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: elements.h:2938
RefineableBuoyantQCrouzeixRaviartElement()
Constructor: call the underlying constructors and initialise the pointer to the Rayleigh number to ad...
unsigned ncont_interpolated_values() const
The total number of continously interpolated values is DIM+1 (DIM fluid velocities and one temperatur...
void disable_ALE()
Final override for disable ALE.
void output(std::ostream &outfile)
Overload the standard output function with the broken default.
virtual void output(std::ostream &outfile)
Output the element data — typically the values at the nodes in a format suitable for post-processing...
Definition: elements.h:2872
void output(std::ostream &outfile)
Overload the standard output function with the broken default.
unsigned num_Z2_flux_terms()
Number of 'flux' terms for Z2 error estimation.
void output(FILE *file_pt)
C-style output function: Broken default.
void further_build()
Call the underlying single-physics element's further_build() functions and make sure that the pointer...
unsigned num_Z2_flux_terms()
The number of Z2 flux terms is the same as that in the fluid element plus that in the advection-diffu...
cstr elem_len * i
Definition: cfortran.h:607
void compute_error(std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
Validate against exact solution at given time Solution is provided via function pointer. Plot at a given number of plot points and compute L2 error and L2 norm of velocity solution over element Overload to broken default.
double dshape_and_dtest_eulerian_at_knot_nst(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Velocity shape and test functions and their derivs w.r.t. to global coords at ipt-th integation point...
virtual void write_tecplot_zone_footer(std::ostream &outfile, const unsigned &nplot) const
Add tecplot zone "footer" to output stream (when plotting nplot points in each "coordinate direction"...
Definition: elements.h:2962
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Fill in the constituent elements' contribution to the residual vector.
void rebuild_from_sons(Mesh *&mesh_pt)
Rebuild from sons: Reconstruct pressure from the (merged) sons This must be specialised for each dime...
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector (wrapper)
void further_build()
Further build: Copy source function pointer from father element.
static double Default_Physical_Constant_Value
The static default value of the Rayleigh number.
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function: Broken default.
void output_fct(std::ostream &outfile, const unsigned &Nplot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
Output function for a time-dependent exact solution. Broken default.
unsigned u_index_adv_diff() const
Overload the index at which the temperature variable is stored. We choose to store is after the fluid...
void get_body_force_nst(const double &time, const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, Vector< double > &result)
Overload the body force in the navier-stokes equations This provides the coupling from the advection-...
virtual unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction") ...
Definition: elements.h:2976
void unfix_pressure(const unsigned &p_dof)
Unpin p_dof-th pressure dof.
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
Definition: elements.h:2458
QAdvectionDiffusionElement elements are linear/quadrilateral/brick-shaped Advection Diffusion element...
const double & pe() const
Peclet number.
void enable_ALE()
(Re-)enable ALE, i.e. take possible mesh motion into account when evaluating the time-derivative. Note: By default, ALE is enabled, at the expense of possibly creating unnecessary work in problems where the mesh is, in fact, stationary.
void output_fct(std::ostream &outfile, const unsigned &Nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output function for an exact solution: Broken default.
char t
Definition: cfortran.h:572
void disable_ALE()
Disable ALE, i.e. assert the mesh is not moving – you do this at your own risk!
virtual double interpolated_p_nst(const Vector< double > &s) const
Return FE interpolated pressure at local coordinate s.
void enable_ALE()
Final override for enable ALE.
double * value_pt(const unsigned &i) const
Return the pointer to the i-the stored value. Typically this is required when direct access to the st...
Definition: nodes.h:322
void enable_ALE()
(Re-)enable ALE, i.e. take possible mesh motion into account when evaluating the time-derivative. Note: By default, ALE is enabled, at the expense of possibly creating unnecessary work in problems where the mesh is, in fact, stationary.
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Compute the element's residual vector and the Jacobian matrix. Jacobian is computed by finite-differe...
void get_interpolated_values(const Vector< double > &s, Vector< double > &values)
Get the continuously interpolated values at the local coordinate s. We choose to put the fluid veloci...
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
HangInfo *const & hanging_pt() const
Return pointer to hanging node data (this refers to the geometric hanging node status) (const version...
Definition: nodes.h:1148
const Vector< double > & g() const
Vector of gravitational components.
void get_wind_adv_diff(const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, Vector< double > &wind) const
Overload the wind function in the advection-diffusion equations. This provides the coupling from the ...
unsigned required_nvalue(const unsigned &n) const
The required number of values stored at the nodes is the sum of the required values of the two single...
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Validate against exact solution. Solution is provided via function pointer. Plot at a given number of...
const double & ra() const
Access function for the Rayleigh number (const version)
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
void fill_in_contribution_to_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
void get_body_force_nst(const double &time, const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, Vector< double > &result)
Overload the body force in the Navier-Stokes equations This provides the coupling from the advection-...
void get_interpolated_values(const unsigned &t, const Vector< double > &s, Vector< double > &values)
Get all continuously interpolated values at the local coordinate s at time level t (t=0: present; t>0...
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element, Call the geometric element's function.
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function: Broken default.
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:733
void disable_ALE()
Final override for disable ALE.
void get_wind_adv_diff(const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, Vector< double > &wind) const
Overload the wind function in the advection-diffusion equations. This provides the coupling from the ...
void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Get 'flux' for Z2 error recovery: Standard flux.from AdvectionDiffusion equations.
void fill_in_contribution_to_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
void fill_in_off_diagonal_jacobian_blocks_analytic(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Compute the contribution of the off-diagonal blocks analytically.
unsigned required_nvalue(const unsigned &n) const
The required number of values stored at the nodes is the sum of the required values of the two single...
void output(FILE *file_pt)
C-style output function: Broken default.
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3841
unsigned nrecovery_order()
The recovery order is that of the NavierStokes elements.
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1207
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. Broken virtual. Needs to be implemented for...
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Validate against exact solution. Solution is provided via function pointer. Plot at a given number of...
double *& ra_pt()
Access function for the pointer to the Rayleigh number.
void further_setup_hanging_nodes()
Perform additional hanging node procedures for variables that are not interpolated by all nodes...
static char t char * s
Definition: cfortran.h:572
std::string scalar_name_paraview(const unsigned &i) const
Name of the i-th scalar field. Default implementation returns V1 for the first one, V2 for the second etc. Can (should!) be overloaded with more meaningful names.
void disable_ALE()
Disable ALE, i.e. assert the mesh is not moving – you do this at your own risk!
void further_setup_hanging_nodes()
The additional hanging node information must be set up for both single-physics elements.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Add the element's contribution to its residual vector and the element Jacobian matrix (wrapper) ...
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:623
void rebuild_from_sons(Mesh *&mesh_pt)
Call the rebuild_from_sons functions for each of the constituent multi-physics elements.
void enable_ALE()
Final override for enable ALE.
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1900
void further_setup_hanging_nodes()
Perform additional hanging node procedures for variables that are not interpolated by all nodes...
Refineable version of Crouzeix Raviart elements. Generic class definitions.
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
double * Ra_pt
Pointer to a new physical variable, the Rayleigh number.
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...
void compute_error(std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
Validate against exact solution at given time Solution is provided via function pointer. Plot at a given number of plot points and compute L2 error and L2 norm of velocity solution over element Call the broken default.
virtual void output_fct(std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output an exact solution over the element.
Definition: elements.h:2911
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Definition: nodes.h:753
double interpolated_u_adv_diff(const Vector< double > &s) const
Return FE representation of function value u(s) at local coordinate s.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
void fill_in_off_diagonal_jacobian_blocks_analytic(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Helper function to get the off-diagonal blocks of the Jacobian matrix analytically.
double Default_Physical_Constant_Value
Default value for physical constants.
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Return the local equation number corresponding to the i-th value at the n-th local node...
Definition: elements.h:1383
unsigned nscalar_paraview() const
Number of scalars/fields output by this element. Broken virtual. Needs to be implemented for each new...
Class that contains data for hanging nodes.
Definition: nodes.h:684
virtual RefineableElement * father_element_pt() const
Return a pointer to the father element.
double *& ra_pt()
Access function for the pointer to the Rayleigh number.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Compute the element's residual Vector.
void output(std::ostream &outfile, const unsigned &nplot)
Output function: x,y,u or x,y,z,u at Nplot^DIM plot points.
BuoyantQCrouzeixRaviartElement()
Constructor: call the underlying constructors and initialise the pointer to the Rayleigh number to po...
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. Broken virtual. Needs to be implemented for...
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 nscalar_paraview() const
Number of scalars/fields output by this element. Broken virtual. Needs to be implemented for each new...
const double & ra() const
Access function for the Rayleigh number (const version)
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Compute the element's residual Vector and the jacobian matrix Virtual function can be overloaded by h...
void output(std::ostream &outfile, const unsigned &nplot)
Output function: Output x, y, u, v, p, theta at Nplot^DIM plot points.
double raw_nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n but do NOT take hanging nodes into account.
Definition: elements.h:2446
unsigned ndof() const
Return the number of equations/dofs in the element.
Definition: elements.h:834
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
void fill_in_contribution_to_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
void rebuild_from_sons(Mesh *&mesh_pt)
Rebuild from sons: empty.
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
void(* UnsteadyExactSolutionFctPt)(const double &, const Vector< double > &, Vector< double > &)
Function pointer for function that computes Vector-valued time-dependent function as ...
Definition: elements.h:1726
virtual std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction") ...
Definition: elements.h:2949
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
Definition: nodes.h:736
virtual void fill_in_contribution_to_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
Add the elemental contribution to the jacobian matrix, mass matrix and the residuals vector...
Definition: elements.cc:1281
void get_Z2_compound_flux_indices(Vector< unsigned > &flux_index)
Fill in which flux components are associated with the fluid measure and which are associated with the...
static double Default_Physical_Constant_Value
The static default value of the Rayleigh number.
unsigned nvertex_node() const
Number of vertex nodes in the element is obtained from the geometric element.
std::string scalar_name_paraview(const unsigned &i) const
Name of the i-th scalar field. Default implementation returns V1 for the first one, V2 for the second etc. Can (should!) be overloaded with more meaningful names.
Refineable version of QAdvectionDiffusionElement. Inherit from the standard QAdvectionDiffusionElemen...
void fill_in_contribution_to_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
Add the element's contribution to its residuals vector, jacobian matrix and mass matrix.
unsigned ncompound_fluxes()
The number of compound fluxes is two (one for the fluid and one for the temperature) ...
unsigned u_index_adv_diff() const
Overload the index at which the temperature variable is stored. We choose to store it after the fluid...
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Calculate the element's contribution to the residual vector. Recall that fill_in_* functions MUST NOT...
void fill_in_off_diagonal_jacobian_blocks_by_fd(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Helper function to get the off-diagonal blocks of the Jacobian matrix by finite differences.
virtual void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Plot the error when compared against a given exact solution . Also calculates the norm of the error a...
Definition: elements.h:2992
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Compute the element's residual Vector and the jacobian matrix using full finite differences, the default implementation.
void output_fct(std::ostream &outfile, const unsigned &Nplot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
Output function for a time-dependent exact solution: Broken default.
static double Default_fd_jacobian_step
Double used for the default finite difference step in elemental jacobian calculations.
Definition: elements.h:1164
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Add the elemental contribution to the jacobian matrix. and the residuals vector. Note that this funct...
Definition: elements.h:1694
A general mesh class.
Definition: mesh.h:74
void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Get the Z2 flux by concatenating the fluxes from the fluid and the advection diffusion elements...
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 * Ra_pt
Pointer to a private data member, the Rayleigh number.
void interpolated_u_nst(const Vector< double > &s, Vector< double > &veloc) const
Compute vector of FE interpolated velocity u at local coordinate s.
virtual unsigned u_index_nst(const unsigned &i) const
Return the index at which the i-th unknown velocity component.