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