multi_domain_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: 1132 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-02-23 05:43:26 +0000 (Tue, 23 Feb 2016) $
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_MULTI_DOMAIN_BOUSSINESQ_ELEMENTS_HEADER
35 #define OOMPH_MULTI_DOMAIN_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 /// Namespace for default parameters in multi-domain Boussinesq
50 //=============================================================
51  namespace MultiDomainBoussinesqHelper
52  {
53 
54  /// Default value for physical constants
56 
57  }
58 /////////////////////////////////////////////////////////////////////////
59 /////////////////////////////////////////////////////////////////////////
60 /////////////////////////////////////////////////////////////////////////
61 
62 
63 //======================nst_bous_class=================================
64 /// Build a refineable Navier Stokes element that inherits from
65 /// ElementWithExternalElement so that it can "communicate" with
66 /// an advection diffusion element that provides the temperature
67 /// in the body force term.
68 //=====================================================================
69 template<class NST_ELEMENT, class AD_ELEMENT>
70 class RefineableNavierStokesBoussinesqElement : public virtual NST_ELEMENT,
71 public virtual ElementWithExternalElement
72 {
73 
74  public:
75 
76  /// \short Constructor: call the underlying constructors and
77  /// initialise the pointer to the Rayleigh number to point
78  /// to the default value of 0.0.
80  NST_ELEMENT(), ElementWithExternalElement()
81  {
83 
84  //There is one interaction: The effect of the advection-diffusion
85  //element onto the buoyancy term
86  this->set_ninteraction(1);
87  }
88 
89  ///Access function for the Rayleigh number (const version)
90  const double &ra() const {return *Ra_pt;}
91 
92  ///Access function for the pointer to the Rayleigh number
93  double* &ra_pt() {return Ra_pt;}
94 
95  /// \short Call the underlying single-physics element's further_build()
96  /// functions and make sure that the pointer to the Rayleigh number
97  /// is passed to the sons. Also make sure that if the external geometric
98  /// Data was ignored in the father it's also ignored in the sons
100  {
101  NST_ELEMENT::further_build();
102 
103  //Cast the pointer to the father element to the specific
104  //element type
106  cast_father_element_pt
108  <NST_ELEMENT,AD_ELEMENT>*>(
109  this->father_element_pt());
110 
111  //Set the pointer to the Rayleigh number to be the same as that in
112  //the father
113  this->Ra_pt = cast_father_element_pt->ra_pt();
114 
115  // Retain ignorance about external geometric data...
116  if (!cast_father_element_pt->external_geometric_data_is_included())
117  {
119  }
120 
121  }
122 
123  /// \short Overload get_body_force_nst() to return the temperature-dependent
124  /// buoyancy force, using the temperature computed by the
125  /// "external" advection diffusion element associated with
126  /// integration point \c ipt.
127  void get_body_force_nst(const double& time, const unsigned& ipt,
128  const Vector<double> &s, const Vector<double> &x,
129  Vector<double> &body_force)
130  {
131  // Set interaction index -- there's only one interaction...
132  const unsigned interaction=0;
133 
134  // Get a pointer to the external element that computes the
135  // the temperature -- we know it's an advection diffusion element.
136  const AD_ELEMENT* adv_diff_el_pt=
137  dynamic_cast<AD_ELEMENT*>(
138  external_element_pt(interaction,ipt));
139 
140  // Get the temperature interpolated from the external element
141  const double interpolated_t =adv_diff_el_pt->
142  interpolated_u_adv_diff(external_element_local_coord(interaction,ipt));
143 
144  // Get vector that indicates the direction of gravity from
145  // the Navier-Stokes equations
146  Vector<double> gravity(NST_ELEMENT::g());
147 
148  // Set the temperature-dependent body force:
149  const unsigned n_dim=this->dim();
150  for (unsigned i=0;i<n_dim;i++)
151  {
152  body_force[i] = -gravity[i]*interpolated_t*ra();
153  }
154 
155  } // end overloaded body force
156 
157 
158  /// \short Compute the element's residual vector and the Jacobian matrix.
160  DenseMatrix<double> &jacobian)
161  {
162 
163  //Get the analytical contribution from the basic Navier-Stokes element
165 
166 #ifdef USE_FD_FOR_DERIVATIVES_WRT_EXTERNAL_DATA_IN_MULTI_DOMAIN_BOUSSINESQ
167 
168  //Get the off-diagonal terms by finite differencing
169  this->fill_in_jacobian_from_external_interaction_by_fd(residuals,jacobian);
170 
171 #else
172 
173  //Get the off-diagonal terms analytically
174  this->fill_in_off_diagonal_block_analytic(residuals,jacobian);
175 
176 #endif
177 
178  }
179 
180 
181 
182  /// Add the element's contribution to its residuals vector,
183  /// jacobian matrix and mass matrix
185  Vector<double> &residuals, DenseMatrix<double> &jacobian,
186  DenseMatrix<double> &mass_matrix)
187  {
188  //Call the standard (Broken) function
189  //which will prevent these elements from being used
190  //in eigenproblems until replaced.
192  residuals,jacobian,mass_matrix);
193  }
194 
195 
196  /// \short Fill in the derivatives of the body force with respect to the
197  /// external unknowns
199  const unsigned& ipt,
200  DenseMatrix<double> &result, Vector<unsigned> &global_eqn_number);
201 
202 
203  /// \short Compute the contribution of the external
204  /// degrees of freedom (temperatures) on the Navier-Stokes equations
206  DenseMatrix<double> &jacobian)
207  {
208  //Local storage for the index in the nodes at which the
209  //Navier-Stokes velocities are stored (we know that this
210  // should be 0,1,2)
211  const unsigned n_dim=this->dim();
212  unsigned u_nodal_nst[n_dim];
213  for(unsigned i=0;i<n_dim;i++)
214  {u_nodal_nst[i] = this->u_index_nst(i);}
215 
216  //Find out how many nodes there are
217  const unsigned n_node = this->nnode();
218 
219  //Set up memory for the shape and test functions and their derivatives
220  Shape psif(n_node), testf(n_node);
221  DShape dpsifdx(n_node,n_dim), dtestfdx(n_node,n_dim);
222 
223  //Number of integration points
224  const unsigned n_intpt = this->integral_pt()->nweight();
225 
226  //Integers to store the local equations and unknowns
227  int local_eqn=0, local_unknown=0;
228 
229  // Local storage for pointers to hang_info objects
230  HangInfo *hang_info_pt=0;
231 
232  //Loop over the integration points
233  for(unsigned ipt=0;ipt<n_intpt;ipt++)
234  {
235  //Get the integral weight
236  double w = this->integral_pt()->weight(ipt);
237 
238  //Call the derivatives of the shape and test functions
239  double J =
240  this->dshape_and_dtest_eulerian_at_knot_nst(ipt,psif,dpsifdx,
241  testf,dtestfdx);
242 
243  //Premultiply the weights and the Jacobian
244  double W = w*J;
245 
246  //Assemble the jacobian terms
247 
248  //Get the derivatives of the body force wrt the unknowns
249  //of the external element
250  DenseMatrix<double> dbody_dexternal_element_data;
251 
252  //Vector of global equation number corresponding to the external
253  //element's data
254  Vector<unsigned> global_eqn_number_of_external_element_data;
255 
256  //Get the appropriate derivatives
258  ipt,dbody_dexternal_element_data,
259  global_eqn_number_of_external_element_data);
260 
261  //Find out how many external data there are
262  const unsigned n_external_element_data =
263  global_eqn_number_of_external_element_data.size();
264 
265  //Loop over the test functions
266  for(unsigned l=0;l<n_node;l++)
267  {
268  //Assemble the contributions of the temperature to
269  //the Navier--Stokes equations (which arise through the buoyancy
270  //body-force term)
271  unsigned n_master = 1;
272  double hang_weight = 1.0;
273 
274  //Local bool (is the node hanging)
275  bool is_node_hanging = this->node_pt(l)->is_hanging();
276 
277  //If the node is hanging, get the number of master nodes
278  if(is_node_hanging)
279  {
280  hang_info_pt = this->node_pt(l)->hanging_pt();
281  n_master = hang_info_pt->nmaster();
282  }
283  //Otherwise there is just one master node, the node itself
284  else
285  {
286  n_master = 1;
287  }
288 
289  //Loop over the master nodes
290  for(unsigned m=0;m<n_master;m++)
291  {
292  //If the node is hanging get weight from master node
293  if(is_node_hanging)
294  {
295  //Get the hang weight from the master node
296  hang_weight = hang_info_pt->master_weight(m);
297  }
298  else
299  {
300  // Node contributes with full weight
301  hang_weight = 1.0;
302  }
303 
304 
305  //Loop over the velocity components in the Navier--Stokes equtions
306  for(unsigned i=0;i<n_dim;i++)
307  {
308  //Get the equation number
309  if(is_node_hanging)
310  {
311  //Get the equation number from the master node
312  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
313  u_nodal_nst[i]);
314  }
315  else
316  {
317  // Local equation number
318  local_eqn = this->nodal_local_eqn(l,u_nodal_nst[i]);
319  }
320 
321  if(local_eqn >= 0)
322  {
323  //Loop over the external data
324  for(unsigned l2=0;l2<n_external_element_data;l2++)
325  {
326  //Find the local equation number corresponding to the global
327  //unknown
328  local_unknown =
329  this->local_eqn_number(
330  global_eqn_number_of_external_element_data[l2]);
331  if(local_unknown >= 0)
332  {
333  //Add contribution to jacobian matrix
334  jacobian(local_eqn,local_unknown)
335  += dbody_dexternal_element_data(i,l2)*testf(l)*hang_weight*W;
336  }
337  }
338  }
339  }
340  }
341  }
342  }
343  }
344 
345  /// Classify dof numbers as in underlying element
346  void get_dof_numbers_for_unknowns(std::list<std::pair<unsigned long,
347  unsigned> >& dof_lookup_list) const
348  {
349  // Call the underlying function
350  NST_ELEMENT::get_dof_numbers_for_unknowns(dof_lookup_list);
351  }
352 
353  /// Get number of dof types from underlying element
354  unsigned ndof_types() const
355  {
356  return NST_ELEMENT::ndof_types();
357  }
358 
359  private:
360 
361  /// Pointer to a private data member, the Rayleigh number
362  double* Ra_pt;
363 
364 };
365 
366 
367 
368 
369 ///////////////////////////////////////////////////////////////////
370 ///////////////////////////////////////////////////////////////////
371 ///////////////////////////////////////////////////////////////////
372 
373 
374 
375 //======================ad_bous_class==================================
376 /// Build an AdvectionDiffusionElement that inherits from
377 /// ElementWithExternalElement so that it can "communicate" with the
378 /// a NavierStokesElement that provides its wind.
379 //=====================================================================
380 template<class AD_ELEMENT, class NST_ELEMENT>
382 public virtual AD_ELEMENT, public virtual ElementWithExternalElement
383 {
384 
385 public:
386 
387  /// \short Constructor: call the underlying constructors
389  AD_ELEMENT(), ElementWithExternalElement()
390  {
391  //There is one interaction
392  this->set_ninteraction(1);
393  }
394 
395 
396  //-----------------------------------------------------------
397  // Note: we're overloading the output functions because the
398  // ===== underlying advection diffusion elements output the
399  // wind which is only available at the integration points
400  // and even then only if the multi domain interaction
401  // has been set up.
402  //-----------------------------------------------------------
403 
404  /// \short Output function:
405  /// Output x, y, theta at Nplot^DIM plot points
406  // Start of output function
407  void output(std::ostream &outfile, const unsigned &nplot)
408  {
409  //vector of local coordinates
410  unsigned n_dim=this->dim();
411  Vector<double> s(n_dim);
412 
413  // Tecplot header info
414  outfile << this->tecplot_zone_string(nplot);
415 
416  // Loop over plot points
417  unsigned num_plot_points=this->nplot_points(nplot);
418  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
419  {
420  // Get local coordinates of plot point
421  this->get_s_plot(iplot,nplot,s);
422 
423  // Output the position of the plot point
424  for(unsigned i=0;i<n_dim;i++)
425  {outfile << this->interpolated_x(s,i) << " ";}
426 
427  // Output the temperature (the advected variable) at the plot point
428  outfile << AD_ELEMENT::interpolated_u_adv_diff(s) << std::endl;
429  }
430  outfile << std::endl;
431 
432  // Write tecplot footer (e.g. FE connectivity lists)
433  this->write_tecplot_zone_footer(outfile,nplot);
434 
435  } //End of output function
436 
437 
438  /// Overload the standard output function with the broken default
439  void output(std::ostream &outfile) {FiniteElement::output(outfile);}
440 
441  /// \short C-style output function: Broken default
442  void output(FILE* file_pt)
443  {FiniteElement::output(file_pt);}
444 
445  /// \short C-style output function: Broken default
446  void output(FILE* file_pt, const unsigned &n_plot)
447  {FiniteElement::output(file_pt,n_plot);}
448 
449 
450  /// \short Overload the wind function in the advection-diffusion equations.
451  /// This provides the coupling from the Navier--Stokes equations to the
452  /// advection-diffusion equations because the wind is the fluid velocity,
453  /// obtained from the "external" element
454  void get_wind_adv_diff(const unsigned& ipt, const Vector<double> &s,
455  const Vector<double>& x, Vector<double>& wind) const
456  {
457  // There is only one interaction
458  unsigned interaction=0;
459 
460  // Dynamic cast "external" element to Navier Stokes element
461  NST_ELEMENT* nst_el_pt= dynamic_cast<NST_ELEMENT*>
462  (external_element_pt(interaction,ipt));
463 
464  //Wind is given by the velocity in the Navier Stokes element
465  nst_el_pt->interpolated_u_nst
466  (external_element_local_coord(interaction,ipt),wind);
467 
468  } //end of get_wind_adv_diff
469 
470 
471  ///\short Compute the element's residual vector and the Jacobian matrix.
473  DenseMatrix<double> &jacobian)
474  {
475  //Get the contribution from the basic advection diffusion element
477 
478 #ifdef USE_FD_FOR_DERIVATIVES_WRT_EXTERNAL_DATA_IN_MULTI_DOMAIN_BOUSSINESQ
479 
480  //Get the off-diagonal terms by finite differencing
481  this->fill_in_jacobian_from_external_interaction_by_fd(residuals,jacobian);
482 
483 #else
484 
485  //Get the off-diagonal terms analytically
486  this->fill_in_off_diagonal_block_analytic(residuals,jacobian);
487 
488 #endif
489 
490  }
491 
492 
493  /// \short Overload the function that must return all field data involved
494  /// in the interaction with the external (Navier Stokes) element.
495  /// Only the velocity dofs in the Navier Stokes element affect the
496  /// interaction with the current element.
498  Vector<std::set<FiniteElement*> > const &external_elements_pt,
499  std::set<std::pair<Data*,unsigned> > &paired_interaction_data);
500 
501  /// \short Add the element's contribution to its residuals vector,
502  /// jacobian matrix and mass matrix
504  Vector<double> &residuals, DenseMatrix<double> &jacobian,
505  DenseMatrix<double> &mass_matrix)
506  {
507  //Call the standard (Broken) function
508  //which will prevent these elements from being used
509  //in eigenproblems until replaced.
511  residuals,jacobian,mass_matrix);
512  }
513 
514 
515  /// \short Fill in the derivatives of the wind with respect to the
516  /// external unknowns
518  const unsigned& ipt, const unsigned &i,
519  Vector<double> &result, Vector<unsigned> &global_eqn_number)
520  {
521  // The interaction index is 0 in this case
522  unsigned interaction=0;
523 
524  // Dynamic cast "other" element to correct type
525  NST_ELEMENT* source_el_pt=dynamic_cast<NST_ELEMENT*>
526  (external_element_pt(interaction,ipt));
527 
528  // Get the external element's derivatives of the velocity with respect
529  // to the data. The wind is just the Navier--Stokes velocity, so this
530  // is all that's required
531  source_el_pt->dinterpolated_u_nst_ddata(
532  external_element_local_coord(interaction,ipt),i,result,
533  global_eqn_number);
534  }
535 
536 
537  /// \short Compute the contribution of the external
538  /// degrees of freedom (velocities) on the advection-diffusion equations
540  DenseMatrix<double> &jacobian)
541  {
542  //Local storage for the index in the nodes at which the temperature
543  //is stored
544  const unsigned u_nodal_adv_diff = this->u_index_adv_diff();
545 
546  //Find out how many nodes there are
547  const unsigned n_node = this->nnode();
548 
549  // Spatial dimension
550  const unsigned n_dim=this->dim();
551 
552  //Set up memory for the shape and test functions and their derivatives
553  Shape psi(n_node), test(n_node);
554  DShape dpsidx(n_node,n_dim), dtestdx(n_node,n_dim);
555 
556  //Number of integration points
557  const unsigned n_intpt = this->integral_pt()->nweight();
558 
559  //Integers to store the local equations and unknowns
560  int local_eqn=0, local_unknown=0;
561 
562  // Local storage for pointers to hang_info objects
563  HangInfo *hang_info_pt=0;
564 
565  //Get the peclet number
566  const double peclet = this->pe();
567 
568  //Loop over the integration points
569  for(unsigned ipt=0;ipt<n_intpt;ipt++)
570  {
571  //Get the integral weight
572  double w = this->integral_pt()->weight(ipt);
573 
574  //Call the derivatives of the shape and test functions
575  double J =
576  this->dshape_and_dtest_eulerian_at_knot_adv_diff(ipt,psi,dpsidx,
577  test,dtestdx);
578 
579  //Premultiply the weights and the Jacobian
580  double W = w*J;
581 
582  //Calculate local values of the derivatives of the solution
583  Vector<double> interpolated_dudx(n_dim,0.0);
584  // Loop over nodes
585  for(unsigned l=0;l<n_node;l++)
586  {
587  // Loop over directions
588  for(unsigned j=0;j<n_dim;j++)
589  {
590  interpolated_dudx[j] +=
591  this->nodal_value(l,u_nodal_adv_diff)*dpsidx(l,j);
592  }
593  }
594 
595  //Get the derivatives of the wind wrt the unknowns
596  //of the external element
597  Vector<double> dwind_dexternal_element_data;
598 
599  //Vector of global equation number corresponding to the external
600  //element's data
601  Vector<unsigned> global_eqn_number_of_external_element_data;
602 
603  //Loop over the wind directions
604  for(unsigned i2=0;i2<n_dim;i2++)
605  {
606  //Get the appropriate derivatives
608  ipt,i2,dwind_dexternal_element_data,
609  global_eqn_number_of_external_element_data);
610 
611 
612  //Find out how many external data there are
613  const unsigned n_external_element_data =
614  global_eqn_number_of_external_element_data.size();
615 
616  //Loop over the test functions
617  for(unsigned l=0;l<n_node;l++)
618  {
619  //Assemble the contributions of the velocities to
620  //the advection-diffusion equations
621  unsigned n_master = 1;
622  double hang_weight = 1.0;
623 
624  //Local bool (is the node hanging)
625  bool is_node_hanging = this->node_pt(l)->is_hanging();
626 
627  //If the node is hanging, get the number of master nodes
628  if(is_node_hanging)
629  {
630  hang_info_pt = this->node_pt(l)->hanging_pt();
631  n_master = hang_info_pt->nmaster();
632  }
633  //Otherwise there is just one master node, the node itself
634  else
635  {
636  n_master = 1;
637  }
638 
639  //Loop over the master nodes
640  for(unsigned m=0;m<n_master;m++)
641  {
642  //If the node is hanging get weight from master node
643  if(is_node_hanging)
644  {
645  //Get the hang weight from the master node
646  hang_weight = hang_info_pt->master_weight(m);
647  }
648  else
649  {
650  // Node contributes with full weight
651  hang_weight = 1.0;
652  }
653 
654  //Get the equation number
655  if(is_node_hanging)
656  {
657  //Get the equation number from the master node
658  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
659  u_nodal_adv_diff);
660  }
661  else
662  {
663  // Local equation number
664  local_eqn = this->nodal_local_eqn(l,u_nodal_adv_diff);
665  }
666 
667  if(local_eqn >= 0)
668  {
669  //Loop over the external data
670  for(unsigned l2=0;l2<n_external_element_data;l2++)
671  {
672  //Find the local equation number corresponding to the global
673  //unknown
674  local_unknown =
675  this->local_eqn_number(
676  global_eqn_number_of_external_element_data[l2]);
677  if(local_unknown >= 0)
678  {
679  //Add contribution to jacobian matrix
680  jacobian(local_eqn,local_unknown)
681  -= peclet*dwind_dexternal_element_data[l2]*
682  interpolated_dudx[i2]*test(l)*hang_weight*W;
683  }
684  }
685  }
686  }
687  }
688  }
689  }
690  }
691 
692  /// Classify dofs for use in block preconditioner
693  void get_dof_numbers_for_unknowns(std::list<std::pair<unsigned long,
694  unsigned> >& dof_lookup_list) const
695  {
696  // number of nodes
697  unsigned n_node = this->nnode();
698 
699  // temporary pair (used to store dof lookup prior to being added to list)
700  std::pair<unsigned,unsigned> dof_lookup;
701 
702  // loop over the nodes
703  for (unsigned n = 0; n < n_node; n++)
704  {
705  // find the number of values at this node
706  unsigned nv = this->node_pt(n)->nvalue();
707 
708  //loop over these values
709  for (unsigned v = 0; v < nv; v++)
710  {
711  // determine local eqn number
712  int local_eqn_number = this->nodal_local_eqn(n, v);
713 
714  // ignore pinned values - far away degrees of freedom resulting
715  // from hanging nodes can be ignored since these are be dealt
716  // with by the element containing their master nodes
717  if (local_eqn_number >= 0)
718  {
719  // store dof lookup in temporary pair: Global equation number
720  // is the first entry in pair
721  dof_lookup.first = this->eqn_number(local_eqn_number);
722 
723  // set dof numbers: Dof number is the second entry in pair
724  dof_lookup.second = 0;
725 
726  // add to list
727  dof_lookup_list.push_front(dof_lookup);
728  }
729  }
730  }
731  }
732 
733 
734  /// Specify number of dof types for use in block preconditioner
735  unsigned ndof_types() const
736  {
737  return 1;
738  }
739 
740 };
741 
742 
743 
744 ///////////////////////////////////////////////////////////////////////////
745 ///////////////////////////////////////////////////////////////////////////
746 ///////////////////////////////////////////////////////////////////////////
747 
748 
749 //================optimised_identification_of_field_data==================
750 /// Overload the function that must return all field data involved
751 /// in the interaction with the external (Navier Stokes) element.
752 /// Only the velocity dofs in the Navier Stokes element affect the
753 /// interaction with the current element.
754 //=======================================================================
755 template<class AD_ELEMENT, class NST_ELEMENT>
758  Vector<std::set<FiniteElement*> > const &external_elements_pt,
759  std::set<std::pair<Data*,unsigned> > &paired_interaction_data)
760  {
761  //There's only one interaction
762  const unsigned interaction = 0;
763 
764  // Loop over each Navier Stokes element in the set of external elements that
765  // affect the current element
766  for(std::set<FiniteElement*>::iterator it=
767  external_elements_pt[interaction].begin();
768  it != external_elements_pt[interaction].end(); it++)
769  {
770 
771  //Cast the external element to a fluid element
772  NST_ELEMENT* external_fluid_el_pt =
773  dynamic_cast<NST_ELEMENT*>(*it);
774 
775  // Loop over the nodes
776  unsigned nnod=external_fluid_el_pt->nnode();
777  for (unsigned j=0;j<nnod;j++)
778  {
779  // Pointer to node (in its incarnation as Data)
780  Data* veloc_data_pt=external_fluid_el_pt->node_pt(j);
781 
782  // Get all velocity dofs
783  const unsigned n_dim=this->dim();
784  for (unsigned i=0;i<n_dim;i++)
785  {
786  // Which value corresponds to the i-th velocity?
787  unsigned val=external_fluid_el_pt->u_index_nst(i);
788 
789  // Turn pointer to Data and index of value into pair
790  // and add to the set
791  paired_interaction_data.insert(std::make_pair(veloc_data_pt,val));
792  }
793  }
794  }
795 } // done
796 
797 
798 /////////////////////////////////////////////////////////////////////////
799 /////////////////////////////////////////////////////////////////////////
800 /////////////////////////////////////////////////////////////////////////
801 
802 
803 //======================class definitions==============================
804 /// Build NavierStokesBoussinesqElement that inherits from
805 /// ElementWithExternalElement so that it can "communicate"
806 /// with AdvectionDiffusionElementWithExternalElement
807 //=====================================================================
808 template<class NST_ELEMENT, class AD_ELEMENT>
810  public virtual NST_ELEMENT, public virtual ElementWithExternalElement
811 {
812 
813 private:
814 
815  /// Pointer to a private data member, the Rayleigh number
816  double* Ra_pt;
817 
818 public:
819 
820  /// \short Constructor: call the underlying constructors and
821  /// initialise the pointer to the Rayleigh number to point
822  /// to the default value of 0.0.
824  {
826 
827  //There is only one interaction
828  this->set_ninteraction(1);
829  }
830 
831  ///Access function for the Rayleigh number (const version)
832  const double &ra() const {return *Ra_pt;}
833 
834  ///Access function for the pointer to the Rayleigh number
835  double* &ra_pt() {return Ra_pt;}
836 
837  ///\short Overload get_body_force_nst to get the temperature "body force"
838  /// from the "source" AdvectionDiffusion element via current integration point
839  void get_body_force_nst(const double& time, const unsigned& ipt,
840  const Vector<double> &s, const Vector<double> &x,
841  Vector<double> &result);
842 
843  /// \short Fill in the derivatives of the body force with respect to the
844  /// external unknowns
846  const unsigned& ipt,
847  DenseMatrix<double> &result, Vector<unsigned> &global_eqn_number);
848 
849 
850  /// \short Compute the element's residual vector and the Jacobian matrix.
851  /// Jacobian is computed by finite-differencing or analytically
853  DenseMatrix<double> &jacobian)
854  {
855 #ifdef USE_FD_JACOBIAN_NST_IN_MULTI_DOMAIN_BOUSSINESQ
856 
857  // This function computes the Jacobian by finite-differencing
859  fill_in_contribution_to_jacobian(residuals,jacobian);
860 
861 #else
862 
863  //Get the contribution from the basic Navier--Stokes element
865 
866  //Get the off-diagonal terms analytically
867  this->fill_in_off_diagonal_block_analytic(residuals,jacobian);
868 
869 #endif
870  }
871 
872  /// \short Add the element's contribution to its residuals vector,
873  /// jacobian matrix and mass matrix
875  Vector<double> &residuals, DenseMatrix<double> &jacobian,
876  DenseMatrix<double> &mass_matrix)
877  {
878  //Call the standard (Broken) function
879  //which will prevent these elements from being used
880  //in eigenproblems until replaced.
882  residuals,jacobian,mass_matrix);
883  }
884 
885  /// \short Compute the contribution of the external
886  /// degrees of freedom (temperatures) on the Navier-Stokes equations
888  DenseMatrix<double> &jacobian)
889  {
890  //Local storage for the index in the nodes at which the
891  //Navier-Stokes velocities are stored (we know that this should be 0,1,2)
892  const unsigned n_dim=this->dim();
893  Vector<unsigned> u_nodal_nst(n_dim);
894  for(unsigned i=0;i<n_dim;i++)
895  {u_nodal_nst[i] = this->u_index_nst(i);}
896 
897  //Find out how many nodes there are
898  const unsigned n_node = this->nnode();
899 
900  //Set up memory for the shape and test functions and their derivatives
901  Shape psif(n_node), testf(n_node);
902  DShape dpsifdx(n_node,n_dim), dtestfdx(n_node,n_dim);
903 
904  //Number of integration points
905  const unsigned n_intpt = this->integral_pt()->nweight();
906 
907  //Integers to store the local equations and unknowns
908  int local_eqn=0, local_unknown=0;
909 
910  //Loop over the integration points
911  for(unsigned ipt=0;ipt<n_intpt;ipt++)
912  {
913  //Get the integral weight
914  double w = this->integral_pt()->weight(ipt);
915 
916  //Call the derivatives of the shape and test functions
917  double J =
918  this->dshape_and_dtest_eulerian_at_knot_nst(ipt,psif,dpsifdx,
919  testf,dtestfdx);
920 
921  //Premultiply the weights and the Jacobian
922  double W = w*J;
923 
924  //Assemble the jacobian terms
925 
926  //Get the derivatives of the body force wrt the unknowns
927  //of the external element
928  DenseMatrix<double> dbody_dexternal_element_data;
929 
930  //Vector of global equation number corresponding to the external
931  //element's data
932  Vector<unsigned> global_eqn_number_of_external_element_data;
933 
934  //Get the appropriate derivatives
936  ipt,dbody_dexternal_element_data,
937  global_eqn_number_of_external_element_data);
938 
939  //Find out how many external data there are
940  const unsigned n_external_element_data =
941  global_eqn_number_of_external_element_data.size();
942 
943  //Loop over the test functions
944  for(unsigned l=0;l<n_node;l++)
945  {
946  //Assemble the contributions of the temperature to
947  //the Navier--Stokes equations (which arise through the buoyancy
948  //body-force term)
949 
950  //Loop over the velocity components in the Navier--Stokes equtions
951  for(unsigned i=0;i<n_dim;i++)
952  {
953  //If it's not a boundary condition
954  local_eqn = this->nodal_local_eqn(l,u_nodal_nst[i]);
955  if(local_eqn >= 0)
956  {
957  //Loop over the external data
958  for(unsigned l2=0;l2<n_external_element_data;l2++)
959  {
960  //Find the local equation number corresponding to the global
961  //unknown
962  local_unknown =
963  this->local_eqn_number(
964  global_eqn_number_of_external_element_data[l2]);
965  if(local_unknown >= 0)
966  {
967  //Add contribution to jacobian matrix
968  jacobian(local_eqn,local_unknown)
969  += dbody_dexternal_element_data(i,l2)*testf(l)*W;
970  }
971  }
972  }
973  }
974  }
975  }
976  }
977 
978 };
979 
980 
981 //============================================================
982 /// Overload get_body_force_nst to get the temperature "body force"
983 /// from the "source" AdvectionDiffusion element via current integration point
984 //========================================================
985 template<class NST_ELEMENT, class AD_ELEMENT>
987 get_body_force_nst(const double& time,
988  const unsigned& ipt,
989  const Vector<double> &s,
990  const Vector<double> &x,
991  Vector<double> &result)
992 {
993  // The interaction index is 0 in this case
994  unsigned interaction=0;
995 
996  // Get the temperature interpolated from the external element
997  const double interpolated_t =
998  dynamic_cast<AD_ELEMENT*>(
999  external_element_pt(interaction,ipt))->
1000  interpolated_u_adv_diff(external_element_local_coord(interaction,ipt));
1001 
1002  // Get vector that indicates the direction of gravity from
1003  // the Navier-Stokes equations
1004  Vector<double> gravity(NST_ELEMENT::g());
1005 
1006  // Temperature-dependent body force:
1007  const unsigned n_dim=this->dim();
1008  for (unsigned i=0;i<n_dim;i++)
1009  {
1010  result[i] = -gravity[i]*interpolated_t*ra();
1011  }
1012 }
1013 
1014 
1015 ///Explicit definition of the face geometry of these elements
1016 template<class NST_ELEMENT, class AD_ELEMENT>
1017  class FaceGeometry<NavierStokesBoussinesqElement<NST_ELEMENT, AD_ELEMENT> > :
1018  public virtual FaceGeometry<NST_ELEMENT>
1019 {
1020  public:
1021  /// \short Constructor calls the constructor of the NST_ELEMENT
1022  /// (Only the Intel compiler seems to need this!)
1023  FaceGeometry() : FaceGeometry<NST_ELEMENT>() {}
1024 
1025  protected:
1026 };
1027 
1028 ///Explicit definition of the face geometry of these elements
1029 template<class NST_ELEMENT, class AD_ELEMENT>
1031  NavierStokesBoussinesqElement<NST_ELEMENT, AD_ELEMENT> > > :
1032  public virtual FaceGeometry<FaceGeometry<NST_ELEMENT> >
1033 {
1034  public:
1035  /// \short Constructor calls the constructor of the NST_ELEMENT
1036  /// (Only the Intel compiler seems to need this!)
1037  FaceGeometry() : FaceGeometry<FaceGeometry<NST_ELEMENT> >() {}
1038 
1039  protected:
1040 };
1041 
1042 
1043 //////////////////////////////////////////////////////////////////////////
1044 //////////////////////////////////////////////////////////////////////////
1045 //////////////////////////////////////////////////////////////////////////
1046 
1047 
1048 //======================class definitions==============================
1049 /// Build AdvectionDiffusionBoussinesqElement that inherits from
1050 /// ElementWithExternalElement so that it can "communicate" with the
1051 /// Navier Stokes element
1052 //=====================================================================
1053 template<class AD_ELEMENT, class NST_ELEMENT>
1055  public virtual AD_ELEMENT, public virtual ElementWithExternalElement
1056 {
1057 
1058 public:
1059 
1060  /// \short Constructor: call the underlying constructors
1063  {
1064  //There is only one interaction
1065  this->set_ninteraction(1);
1066  }
1067 
1068  /// \short Overload the wind function in the advection-diffusion equations.
1069  /// This provides the coupling from the Navier--Stokes equations to the
1070  /// advection-diffusion equations because the wind is the fluid velocity,
1071  /// obtained from the source element in the other mesh
1072  void get_wind_adv_diff(const unsigned& ipt, const Vector<double> &s,
1073  const Vector<double>& x, Vector<double>& wind) const;
1074 
1075 
1076  //-----------------------------------------------------------
1077  // Note: we're overloading the output functions because the
1078  // ===== underlying advection diffusion elements output the
1079  // wind which is only available at the integration points
1080  // and even then only if the multi domain interaction
1081  // has been set up.
1082  //-----------------------------------------------------------
1083 
1084  /// \short Output function:
1085  /// Output x, y, theta at Nplot^DIM plot points
1086  // Start of output function
1087  void output(std::ostream &outfile, const unsigned &nplot)
1088  {
1089  // Element dimension
1090  const unsigned n_dim=this->dim();
1091 
1092  //vector of local coordinates
1093  Vector<double> s(n_dim);
1094 
1095  // Tecplot header info
1096  outfile << this->tecplot_zone_string(nplot);
1097 
1098  // Loop over plot points
1099  unsigned num_plot_points=this->nplot_points(nplot);
1100  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
1101  {
1102  // Get local coordinates of plot point
1103  this->get_s_plot(iplot,nplot,s);
1104 
1105  // Output the position of the plot point
1106  for(unsigned i=0;i<n_dim;i++)
1107  {outfile << this->interpolated_x(s,i) << " ";}
1108 
1109  // Output the temperature (the advected variable) at the plot point
1110  outfile << AD_ELEMENT::interpolated_u_adv_diff(s) << std::endl;
1111  }
1112  outfile << std::endl;
1113 
1114  // Write tecplot footer (e.g. FE connectivity lists)
1115  this->write_tecplot_zone_footer(outfile,nplot);
1116 
1117  } //End of output function
1118 
1119 
1120  /// Overload the standard output function with the broken default
1121  void output(std::ostream &outfile) {FiniteElement::output(outfile);}
1122 
1123  /// \short C-style output function: Broken default
1124  void output(FILE* file_pt)
1125  {FiniteElement::output(file_pt);}
1126 
1127  /// \short C-style output function: Broken default
1128  void output(FILE* file_pt, const unsigned &n_plot)
1129  {FiniteElement::output(file_pt,n_plot);}
1130 
1131 
1132 
1133  /// \short Fill in the derivatives of the wind with respect to the
1134  /// external unknowns
1136  const unsigned& ipt, const unsigned &i,
1137  Vector<double> &result, Vector<unsigned> &global_eqn_number);
1138 
1139 
1140  ///\short Compute the element's residual vector and the Jacobian matrix.
1141  /// Jacobian is computed by finite-differencing.
1143  DenseMatrix<double> &jacobian)
1144  {
1145 #ifdef USE_FD_JACOBIAN_IN_MULTI_DOMAIN_BOUSSINESQ
1146 
1147  // This function computes the Jacobian by finite-differencing
1149  fill_in_contribution_to_jacobian(residuals,jacobian);
1150 
1151 #else
1152 
1153  //Get the contribution from the basic advection-diffusion element
1155 
1156  //Get the off-diagonal terms analytically
1157  this->fill_in_off_diagonal_block_analytic(residuals,jacobian);
1158 
1159 #endif
1160  }
1161 
1162  /// \short Add the element's contribution to its residuals vector,
1163  /// jacobian matrix and mass matrix
1165  Vector<double> &residuals, DenseMatrix<double> &jacobian,
1166  DenseMatrix<double> &mass_matrix)
1167  {
1168  //Call the standard (Broken) function
1169  //which will prevent these elements from being used
1170  //in eigenproblems until replaced.
1172  residuals,jacobian,mass_matrix);
1173  }
1174 
1175 
1176  /// \short Compute the contribution of the external
1177  /// degrees of freedom (velocities) on the AdvectionDiffusion equations
1179  DenseMatrix<double> &jacobian)
1180  {
1181  //Local storage for the index at which the temperature is stored
1182  const unsigned u_nodal_adv_diff = this->u_index_adv_diff();
1183 
1184  //Find out how many nodes there are
1185  const unsigned n_node = this->nnode();
1186 
1187  // Element dimension
1188  const unsigned n_dim=this->dim();
1189 
1190  //Set up memory for the shape and test functions and their derivatives
1191  Shape psi(n_node), test(n_node);
1192  DShape dpsidx(n_node,n_dim), dtestdx(n_node,n_dim);
1193 
1194  //Number of integration points
1195  const unsigned n_intpt = this->integral_pt()->nweight();
1196 
1197  //Integers to store the local equations and unknowns
1198  int local_eqn=0, local_unknown=0;
1199 
1200  //Get the peclet number
1201  const double peclet = this->pe();
1202 
1203  //Loop over the integration points
1204  for(unsigned ipt=0;ipt<n_intpt;ipt++)
1205  {
1206  //Get the integral weight
1207  double w = this->integral_pt()->weight(ipt);
1208 
1209  //Call the derivatives of the shape and test functions
1210  double J =
1211  this->dshape_and_dtest_eulerian_at_knot_adv_diff(ipt,psi,dpsidx,
1212  test,dtestdx);
1213 
1214  //Premultiply the weights and the Jacobian
1215  double W = w*J;
1216 
1217  //Calculate local values of the derivatives of the solution
1218  Vector<double> interpolated_dudx(n_dim,0.0);
1219  // Loop over nodes
1220  for(unsigned l=0;l<n_node;l++)
1221  {
1222  // Loop over directions
1223  for(unsigned j=0;j<n_dim;j++)
1224  {
1225  interpolated_dudx[j] +=
1226  this->raw_nodal_value(l,u_nodal_adv_diff)*dpsidx(l,j);
1227  }
1228  }
1229 
1230  //Get the derivatives of the wind wrt the unknowns
1231  //of the external element
1232  Vector<double> dwind_dexternal_element_data;
1233 
1234  //Vector of global equation number corresponding to the external
1235  //element's data
1236  Vector<unsigned> global_eqn_number_of_external_element_data;
1237 
1238 
1239  //Loop over the wind directions
1240  for(unsigned i2=0;i2<n_dim;i2++)
1241  {
1242  //Get the appropriate derivatives
1244  ipt,i2,dwind_dexternal_element_data,
1245  global_eqn_number_of_external_element_data);
1246 
1247  //Find out how many external data there are
1248  const unsigned n_external_element_data =
1249  global_eqn_number_of_external_element_data.size();
1250 
1251  //Loop over the test functions
1252  for(unsigned l=0;l<n_node;l++)
1253  {
1254  //Assemble the contributions of the velocities
1255  //the Advection-Diffusion equations
1256 
1257  //If it's not a boundary condition
1258  local_eqn = this->nodal_local_eqn(l,u_nodal_adv_diff);
1259  if(local_eqn >= 0)
1260  {
1261  //Loop over the external data
1262  for(unsigned l2=0;l2<n_external_element_data;l2++)
1263  {
1264  //Find the local equation number corresponding to the global
1265  //unknown
1266  local_unknown =
1267  this->local_eqn_number(
1268  global_eqn_number_of_external_element_data[l2]);
1269  if(local_unknown >= 0)
1270  {
1271  //Add contribution to jacobian matrix
1272  jacobian(local_eqn,local_unknown)
1273  -= peclet*dwind_dexternal_element_data[l2]*
1274  interpolated_dudx[i2]*test(l)*W;
1275  }
1276  }
1277  }
1278  }
1279  }
1280  }
1281  }
1282 
1283 };
1284 
1285 
1286 
1287 //==========================================================================
1288 /// \short Overload the wind function in the advection-diffusion equations.
1289 /// This provides the coupling from the Navier--Stokes equations to the
1290 /// advection-diffusion equations because the wind is the fluid velocity,
1291 /// obtained from the source elements in the other mesh
1292 //==========================================================================
1293 template<class AD_ELEMENT, class NST_ELEMENT>
1296 (const unsigned& ipt,const Vector<double> &s,const Vector<double>& x,
1297  Vector<double>& wind) const
1298 {
1299  // The interaction index is 0 in this case
1300  unsigned interaction=0;
1301 
1302  // Dynamic cast "other" element to correct type
1303  NST_ELEMENT* source_el_pt=
1304  dynamic_cast<NST_ELEMENT*>(external_element_pt(interaction,ipt));
1305 
1306  //The wind function is simply the velocity at the points of the "other" el
1307  source_el_pt->interpolated_u_nst
1308  (external_element_local_coord(interaction,ipt),wind);
1309 }
1310 
1311 //=========================================================================
1312 /// Fill in the derivatives of the wind with respect to the external
1313 /// unknowns in the advection-diffusion equations
1314 //=========================================================================
1315 template<class AD_ELEMENT, class NST_ELEMENT>
1318  const unsigned &i,
1319  Vector<double> &result,
1320  Vector<unsigned> &global_eqn_number)
1321 {
1322  // The interaction index is 0 in this case
1323  unsigned interaction=0;
1324 
1325  // Dynamic cast "other" element to correct type
1326  NST_ELEMENT* source_el_pt=
1327  dynamic_cast<NST_ELEMENT*>(external_element_pt(interaction,ipt));
1328 
1329  // Get the external element's derivatives of the velocity with respect
1330  // to the data. The wind is just the Navier--Stokes velocity, so this
1331  // is all that's required
1332  source_el_pt->dinterpolated_u_nst_ddata(
1333  external_element_local_coord(interaction,ipt),i,result,
1334  global_eqn_number);
1335 }
1336 
1337 //////////////////////////////////////////////////////////////////////
1338 //////////////////////////////////////////////////////////////////////
1339 //////////////////////////////////////////////////////////////////////
1340 
1341 //=========================================================================
1342 /// Fill in the derivatives of the body force with respect to the external
1343 /// unknowns in the Navier--Stokes equations
1344 //=========================================================================
1345 template<class NST_ELEMENT, class AD_ELEMENT>
1348  DenseMatrix<double> &result,
1349  Vector<unsigned> &global_eqn_number)
1350 {
1351  // Get vector that indicates the direction of gravity from
1352  // the Navier-Stokes equations
1353  Vector<double> gravity(NST_ELEMENT::g());
1354 
1355  // The interaction index is 0 in this case
1356  unsigned interaction=0;
1357 
1358  // Get the external element's derivatives
1359  Vector<double> du_adv_diff_ddata;
1360 
1361  // Dynamic cast "other" element to correct type
1362  AD_ELEMENT* source_el_pt=
1363  dynamic_cast<AD_ELEMENT*>(external_element_pt(interaction,ipt));
1364 #ifdef PARANOID
1365  if (source_el_pt==0)
1366  {
1367  throw OomphLibError(
1368  "External element could not be cast to AD_ELEMENT\n",
1369  OOMPH_CURRENT_FUNCTION,
1370  OOMPH_EXCEPTION_LOCATION);
1371  }
1372 #endif
1373 
1374  // Get derivatives
1375  source_el_pt->dinterpolated_u_adv_diff_ddata(
1376  external_element_local_coord(interaction,ipt),du_adv_diff_ddata,
1377  global_eqn_number);
1378 
1379  //Find the number of external data
1380  unsigned n_external_element_data = du_adv_diff_ddata.size();
1381 
1382  //Set the size of the matrix to be returned
1383  const unsigned n_dim=this->dim();
1384  result.resize(n_dim,n_external_element_data);
1385 
1386  // Temperature-dependent body force:
1387  for (unsigned i=0;i<n_dim;i++)
1388  {
1389  //Loop over the external data
1390  for(unsigned n=0;n<n_external_element_data;n++)
1391  {
1392  result(i,n) = -gravity[i]*du_adv_diff_ddata[n]*ra();
1393  }
1394  }
1395 }
1396 
1397 
1398 //==========start_of_get_dbody_force=========== ===========================
1399 /// Fill in the derivatives of the body force with respect to the external
1400 /// unknowns in the Navier--Stokes equations
1401 //=========================================================================
1402 template<class NST_ELEMENT, class AD_ELEMENT>
1405  DenseMatrix<double> &result,
1406  Vector<unsigned> &global_eqn_number)
1407 {
1408  // Get vector that indicates the direction of gravity from
1409  // the Navier-Stokes equations
1410  Vector<double> gravity(NST_ELEMENT::g());
1411 
1412  // The interaction index is 0 in this case
1413  unsigned interaction=0;
1414 
1415  // Get the external element's derivatives
1416  Vector<double> du_adv_diff_ddata;
1417 
1418  // Dynamic cast "other" element to correct type
1419  AD_ELEMENT* source_el_pt=
1420  dynamic_cast<AD_ELEMENT*>(external_element_pt(interaction,ipt));
1421 #ifdef PARANOID
1422  if (source_el_pt==0)
1423  {
1424  throw OomphLibError(
1425  "External element could not be cast to AD_ELEMENT\n",
1426  OOMPH_CURRENT_FUNCTION,
1427  OOMPH_EXCEPTION_LOCATION);
1428  }
1429 #endif
1430 
1431  // Get derivatives
1432  source_el_pt->dinterpolated_u_adv_diff_ddata(
1433  external_element_local_coord(interaction,ipt),du_adv_diff_ddata,
1434  global_eqn_number);
1435 
1436  //Find the number of external data
1437  unsigned n_external_element_data = du_adv_diff_ddata.size();
1438 
1439  //Set the size of the matrix to be returned
1440  const unsigned n_dim=this->dim();
1441  result.resize(n_dim,n_external_element_data);
1442 
1443  // Temperature-dependent body force:
1444  for (unsigned i=0;i<n_dim;i++)
1445  {
1446  //Loop over the external data
1447  for(unsigned n=0;n<n_external_element_data;n++)
1448  {
1449  result(i,n) = -gravity[i]*du_adv_diff_ddata[n]*ra();
1450  }
1451  }
1452 
1453 } // end_of_get_dbody_force
1454 
1455 }
1456 
1457 #endif
void get_dof_numbers_for_unknowns(std::list< std::pair< unsigned long, unsigned > > &dof_lookup_list) const
Classify dofs for use in block preconditioner.
void output(std::ostream &outfile, const unsigned &nplot)
Output function: Output x, y, theta at Nplot^DIM plot points.
double * Ra_pt
Pointer to a private data member, the Rayleigh number.
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...
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
void get_dwind_adv_diff_dexternal_element_data(const unsigned &ipt, const unsigned &i, Vector< double > &result, Vector< unsigned > &global_eqn_number)
Fill in the derivatives of the wind with respect to the external unknowns.
void identify_all_field_data_for_external_interaction(Vector< std::set< FiniteElement * > > const &external_elements_pt, std::set< std::pair< Data *, unsigned > > &paired_interaction_data)
Overload the function that must return all field data involved in the interaction with the external (...
void get_dwind_adv_diff_dexternal_element_data(const unsigned &ipt, const unsigned &i, Vector< double > &result, Vector< unsigned > &global_eqn_number)
Fill in the derivatives of the wind with respect to the external unknowns.
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.
Vector< double > & external_element_local_coord(const unsigned &interaction_index, const unsigned &ipt)
Access function to get source element's local coords for specified interaction index at specified int...
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 fill_in_contribution_to_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
void get_dbody_force_nst_dexternal_element_data(const unsigned &ipt, DenseMatrix< double > &result, Vector< unsigned > &global_eqn_number)
Fill in the derivatives of the body force with respect to the external unknowns.
RefineableNavierStokesBoussinesqElement()
Constructor: call the underlying constructors and initialise the pointer to the Rayleigh number to po...
cstr elem_len * i
Definition: cfortran.h:607
FaceGeometry()
Constructor calls the constructor of the NST_ELEMENT (Only the Intel compiler seems to need this!) ...
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_jacobian_from_external_interaction_by_fd(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Calculate the contributions to the jacobian from all external interaction degrees of freedom (geometr...
int local_eqn_number(const unsigned long &ieqn_global) const
Return the local equation number corresponding to the ieqn_global-th global equation number...
Definition: elements.h:731
bool external_geometric_data_is_included() const
Is the external geometric data taken into account when forming the Jacobian?
const double & ra() const
Access function for the Rayleigh number (const version)
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
unsigned long eqn_number(const unsigned &ieqn_local) const
Return the global equation number corresponding to the ieqn_local-th local equation number...
Definition: elements.h:709
void output(FILE *file_pt)
C-style output function: Broken default.
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
RefineableAdvectionDiffusionBoussinesqElement()
Constructor: call the underlying constructors.
double * Ra_pt
Pointer to a private data member, the Rayleigh number.
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.
void set_ninteraction(const unsigned &n_interaction)
Set the number of interactions in the element This function is usually called in the specific element...
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
void output(std::ostream &outfile)
Overload the standard output function with the broken default.
FaceGeometry()
Constructor calls the constructor of the NST_ELEMENT (Only the Intel compiler seems to need this!) ...
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:448
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function: Broken default.
const double & pe() const
Peclet number.
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:733
void get_dbody_force_nst_dexternal_element_data(const unsigned &ipt, DenseMatrix< double > &result, Vector< unsigned > &global_eqn_number)
Fill in the derivatives of the body force with respect to the external unknowns.
const double & ra() const
Access function for the Rayleigh number (const version)
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 ...
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
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1207
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) ...
double *& ra_pt()
Access function for the pointer to the Rayleigh number.
void further_build()
Call the underlying single-physics element's further_build() functions and make sure that the pointer...
static char t char * s
Definition: cfortran.h:572
NavierStokesBoussinesqElement()
Constructor: call the underlying constructors and initialise the pointer to the Rayleigh number to po...
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1900
unsigned ndof_types() const
Specify number of dof types for use in block preconditioner.
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:89
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Definition: nodes.h:753
AdvectionDiffusionBoussinesqElement()
Constructor: call the underlying constructors.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
double Default_Physical_Constant_Value
Default value for physical constants.
void fill_in_off_diagonal_block_analytic(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Compute the contribution of the external degrees of freedom (velocities) on the advection-diffusion e...
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
Class that contains data for hanging nodes.
Definition: nodes.h:684
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Compute the element's residual vector and the Jacobian matrix.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Compute the element's residual vector and the Jacobian matrix.
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
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 fill_in_off_diagonal_block_analytic(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Compute the contribution of the external degrees of freedom (temperatures) on the Navier-Stokes equat...
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2470
FiniteElement *& external_element_pt(const unsigned &interaction_index, const unsigned &ipt)
Access function to source element for specified interaction index at specified integration point...
void output(FILE *file_pt)
C-style output function: Broken default.
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.
void ignore_external_geometric_data()
Do not include any external geometric data when computing the element's Jacobian. This function shoul...
void get_dof_numbers_for_unknowns(std::list< std::pair< unsigned long, unsigned > > &dof_lookup_list) const
Classify dof numbers as in underlying element.
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
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
void get_body_force_nst(const double &time, const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, Vector< double > &body_force)
Overload get_body_force_nst() to return the temperature-dependent buoyancy force, using the temperatu...
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...
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
double *& ra_pt()
Access function for the pointer to the Rayleigh number.
void fill_in_off_diagonal_block_analytic(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Compute the contribution of the external degrees of freedom (velocities) on the AdvectionDiffusion eq...
void fill_in_off_diagonal_block_analytic(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Compute the contribution of the external degrees of freedom (temperatures) on the Navier-Stokes equat...
void get_body_force_nst(const double &time, const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, Vector< double > &result)
Overload get_body_force_nst to get the temperature "body force" from the "source" AdvectionDiffusion ...
void output(std::ostream &outfile)
Overload the standard output function with the broken default.
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function: Broken default.
unsigned ndof_types() const
Get number of dof types from underlying element.
void output(std::ostream &outfile, const unsigned &nplot)
Output function: Output x, y, theta at Nplot^DIM plot points.
void resize(const unsigned long &n)
Definition: matrices.h:511