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-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26// Header 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
40namespace 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 {
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.
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
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.
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
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.
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
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.
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
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
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
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.
898 residuals, jacobian, mass_matrix);
899 }
900
901 /// Compute the contribution of the external
902 /// degrees of freedom (temperatures) on the Navier-Stokes equations
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>
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!)
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.
1170 DenseMatrix<double>& jacobian)
1171 {
1172#ifdef USE_FD_JACOBIAN_IN_MULTI_DOMAIN_BOUSSINESQ
1173
1174 // This function computes the Jacobian by finite-differencing
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.
1200 residuals, jacobian, mass_matrix);
1201 }
1202
1203
1204 /// Compute the contribution of the external
1205 /// degrees of freedom (velocities) on the AdvectionDiffusion equations
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
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
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.
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:86
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:483
void resize(const unsigned long &n)
Resize to a square nxn matrix; any values already present will be transfered.
Definition: matrices.h:498
This is a base class for all elements that require external sources (e.g. FSI, multi-domain problems ...
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...
void set_ninteraction(const unsigned &n_interaction)
Set the number of interactions in the element This function is usually called in the specific element...
bool external_geometric_data_is_included() const
Is the external geometric data taken into account when forming the Jacobian?
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...
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 ignore_external_geometric_data()
Do not include any external geometric data when computing the element's Jacobian. This function shoul...
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in the element's contribution to the Jacobian matrix and the residual vector: Done by finite dif...
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!...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:4998
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
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:2593
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:3050
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:3161
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:3962
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:1432
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2611
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &shifted_to_interior=false) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: elements.h:3148
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:3186
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
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:2576
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:3174
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:1322
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:704
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:726
Class that contains data for hanging nodes.
Definition: nodes.h:742
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
Definition: nodes.h:791
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:785
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Definition: nodes.h:808
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
////////////////////////////////////////////////////////////////////// //////////////////////////////...
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 ...
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...
double *& ra_pt()
Access function for the pointer to 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 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)
NavierStokesBoussinesqElement()
Constructor: call the underlying constructors and initialise the pointer to the Rayleigh number to po...
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1285
HangInfo *const & hanging_pt() const
Return pointer to hanging node data (this refers to the geometric hanging node status) (const version...
Definition: nodes.h:1228
An OomphLibError object which should be thrown when an run-time error is encountered....
//////////////////////////////////////////////////////////////// ////////////////////////////////////...
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 get_dof_numbers_for_unknowns(std::list< std::pair< unsigned long, unsigned > > &dof_lookup_list) const
Classify dofs for use in block preconditioner.
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_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.
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 get_dof_numbers_for_unknowns(std::list< std::pair< unsigned long, unsigned > > &dof_lookup_list) const
Classify dof numbers as in underlying element.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Compute the element's residual vector and the Jacobian matrix.
double *& ra_pt()
Access function for the pointer to the Rayleigh number.
RefineableNavierStokesBoussinesqElement()
Constructor: call the underlying constructors and initialise the pointer to the Rayleigh number to po...
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
double Default_Physical_Constant_Value
Default value for physical constants.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...