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_BOUSSINESQ_ELEMENTS_HEADER
31#define OOMPH_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 /// //////////////////////////////////////////////////////////////////////
44 /// //////////////////////////////////////////////////////////////////////
45
46
47 //======================class definition==============================
48 /// A class that solves the Boussinesq approximation of the Navier--Stokes
49 /// and energy equations by coupling two pre-existing classes.
50 /// The QAdvectionDiffusionElement with bi-quadratic interpolation for the
51 /// scalar variable (temperature) and
52 /// QCrouzeixRaviartElement which solves the Navier--Stokes equations
53 /// using bi-quadratic interpolation for the velocities and a discontinuous
54 /// bi-linear interpolation for the pressure. Note that we are free to
55 /// choose the order in which we store the variables at the nodes. In this
56 /// case we choose to store the variables in the order fluid velocities
57 /// followed by temperature. We must, therefore, overload the function
58 /// AdvectionDiffusionEquations<DIM>::u_index_adv_diff() to indicate that
59 /// the temperature is stored at the DIM-th position not the 0-th. We do not
60 /// need to overload the corresponding function in the
61 /// NavierStokesEquations<DIM> class because the velocities are stored
62 /// first.
63 //=========================================================================
64 template<unsigned DIM>
66 : public virtual QAdvectionDiffusionElement<DIM, 3>,
67 public virtual QCrouzeixRaviartElement<DIM>
68 {
69 private:
70 /// Pointer to a private data member, the Rayleigh number
71 double* Ra_pt;
72
73 /// The static default value of the Rayleigh number
75
76 public:
77 /// Constructor: call the underlying constructors and
78 /// initialise the pointer to the Rayleigh number to point
79 /// to the default value of 0.0.
81 : QAdvectionDiffusionElement<DIM, 3>(), QCrouzeixRaviartElement<DIM>()
82 {
84 }
85
86 /// Unpin p_dof-th pressure dof
87 void unfix_pressure(const unsigned& p_dof)
88 {
89 this->internal_data_pt(this->P_nst_internal_index)->unpin(p_dof);
90 }
91
92 /// The required number of values stored at the nodes is the sum of
93 /// the required values of the two single-physics elements. Note that this
94 /// step is generic for any multi-physics element of this type.
95 unsigned required_nvalue(const unsigned& n) const
96 {
97 return (QAdvectionDiffusionElement<DIM, 3>::required_nvalue(n) +
98 QCrouzeixRaviartElement<DIM>::required_nvalue(n));
99 }
100
101 /// Access function for the Rayleigh number (const version)
102 const double& ra() const
103 {
104 return *Ra_pt;
105 }
106
107 /// Access function for the pointer to the Rayleigh number
108 double*& ra_pt()
109 {
110 return Ra_pt;
111 }
112
113 /// Final override for disable ALE
115 {
116 // Disable ALE in both sets of equations
117 NavierStokesEquations<DIM>::disable_ALE();
118 AdvectionDiffusionEquations<DIM>::disable_ALE();
119 }
120
121 /// Final override for enable ALE
123 {
124 // Enable ALE in both sets of equations
125 NavierStokesEquations<DIM>::enable_ALE();
126 AdvectionDiffusionEquations<DIM>::enable_ALE();
127 }
128
129 /// Number of scalars/fields output by this element. Broken
130 /// virtual. Needs to be implemented for each new specific element type.
131 /// Temporary dummy
132 unsigned nscalar_paraview() const
133 {
134 throw OomphLibError(
135 "This function hasn't been implemented for this element",
136 OOMPH_CURRENT_FUNCTION,
137 OOMPH_EXCEPTION_LOCATION);
138
139 // Dummy unsigned
140 return 0;
141 }
142
143 /// Write values of the i-th scalar field at the plot points. Broken
144 /// virtual. Needs to be implemented for each new specific element type.
145 /// Temporary dummy
146 void scalar_value_paraview(std::ofstream& file_out,
147 const unsigned& i,
148 const unsigned& nplot) const
149 {
150 throw OomphLibError(
151 "This function hasn't been implemented for this element",
152 OOMPH_CURRENT_FUNCTION,
153 OOMPH_EXCEPTION_LOCATION);
154 }
155
156
157 /// Name of the i-th scalar field. Default implementation
158 /// returns V1 for the first one, V2 for the second etc. Can (should!) be
159 /// overloaded with more meaningful names.
160 std::string scalar_name_paraview(const unsigned& i) const
161 {
162 return "V" + StringConversion::to_string(i);
163 }
164
165
166 /// Overload the standard output function with the broken default
167 void output(std::ostream& outfile)
168 {
169 FiniteElement::output(outfile);
170 }
171
172 /// Output function:
173 /// Output x, y, u, v, p, theta at Nplot^DIM plot points
174 // Start of output function
175 void output(std::ostream& outfile, const unsigned& nplot)
176 {
177 // vector of local coordinates
178 Vector<double> s(DIM);
179
180 // Tecplot header info
181 outfile << this->tecplot_zone_string(nplot);
182
183 // Loop over plot points
184 unsigned num_plot_points = this->nplot_points(nplot);
185 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
186 {
187 // Get local coordinates of plot point
188 this->get_s_plot(iplot, nplot, s);
189
190 // Output the position of the plot point
191 for (unsigned i = 0; i < DIM; i++)
192 {
193 outfile << this->interpolated_x(s, i) << " ";
194 }
195
196 // Output the fluid velocities at the plot point
197 for (unsigned i = 0; i < DIM; i++)
198 {
199 outfile << this->interpolated_u_nst(s, i) << " ";
200 }
201
202 // Output the fluid pressure at the plot point
203 outfile << this->interpolated_p_nst(s) << " ";
204
205 // Output the temperature (the advected variable) at the plot point
206 outfile << this->interpolated_u_adv_diff(s) << std::endl;
207 }
208 outfile << std::endl;
209
210 // Write tecplot footer (e.g. FE connectivity lists)
211 this->write_tecplot_zone_footer(outfile, nplot);
212 } // End of output function
213
214
215 /// C-style output function: Broken default
216 void output(FILE* file_pt)
217 {
218 FiniteElement::output(file_pt);
219 }
220
221 /// C-style output function: Broken default
222 void output(FILE* file_pt, const unsigned& n_plot)
223 {
224 FiniteElement::output(file_pt, n_plot);
225 }
226
227 /// Output function for an exact solution: Broken default
228 void output_fct(std::ostream& outfile,
229 const unsigned& Nplot,
230 FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
231 {
232 FiniteElement::output_fct(outfile, Nplot, exact_soln_pt);
233 }
234
235
236 /// Output function for a time-dependent exact solution:
237 /// Broken default.
238 void output_fct(std::ostream& outfile,
239 const unsigned& Nplot,
240 const double& time,
241 FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
242 {
243 FiniteElement::output_fct(outfile, Nplot, time, exact_soln_pt);
244 }
245
246 /// Overload the index at which the temperature
247 /// variable is stored. We choose to store it after the fluid velocities.
248 inline unsigned u_index_adv_diff() const
249 {
250 return DIM;
251 }
252
253 /// Validate against exact solution at given time
254 /// Solution is provided via function pointer.
255 /// Plot at a given number of plot points and compute L2 error
256 /// and L2 norm of velocity solution over element
257 /// Call the broken default
258 void compute_error(std::ostream& outfile,
259 FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt,
260 const double& time,
261 double& error,
262 double& norm)
263 {
264 FiniteElement::compute_error(outfile, exact_soln_pt, time, error, norm);
265 }
266
267 /// Validate against exact solution.
268 /// Solution is provided via function pointer.
269 /// Plot at a given number of plot points and compute L2 error
270 /// and L2 norm of velocity solution over element
271 /// Call the broken default
272 void compute_error(std::ostream& outfile,
273 FiniteElement::SteadyExactSolutionFctPt exact_soln_pt,
274 double& error,
275 double& norm)
276 {
277 FiniteElement::compute_error(outfile, exact_soln_pt, error, norm);
278 }
279
280 /// Overload the wind function in the advection-diffusion equations.
281 /// This provides the coupling from the Navier--Stokes equations to the
282 /// advection-diffusion equations because the wind is the fluid velocity.
283 void get_wind_adv_diff(const unsigned& ipt,
284 const Vector<double>& s,
285 const Vector<double>& x,
286 Vector<double>& wind) const
287 {
288 // The wind function is simply the velocity at the points
289 this->interpolated_u_nst(s, wind);
290 }
291
292
293 /// Overload the body force in the Navier-Stokes equations
294 /// This provides the coupling from the advection-diffusion equations
295 /// to the Navier--Stokes equations, the body force is the
296 /// temperature multiplied by the Rayleigh number acting in the
297 /// direction opposite to gravity.
298 void get_body_force_nst(const double& time,
299 const unsigned& ipt,
300 const Vector<double>& s,
301 const Vector<double>& x,
302 Vector<double>& result)
303 {
304 // Get the temperature
305 const double interpolated_t = this->interpolated_u_adv_diff(s);
306
307 // Get vector that indicates the direction of gravity from
308 // the Navier-Stokes equations
309 Vector<double> gravity(NavierStokesEquations<DIM>::g());
310
311 // Temperature-dependent body force:
312 for (unsigned i = 0; i < DIM; i++)
313 {
314 result[i] = -gravity[i] * interpolated_t * ra();
315 }
316 } // end of get_body_force
317
318 /// Calculate the element's contribution to the residual vector.
319 /// Recall that fill_in_* functions MUST NOT initialise the entries
320 /// in the vector to zero. This allows us to call the
321 /// fill_in_* functions of the constituent single-physics elements
322 /// sequentially, without wiping out any previously computed entries.
323 void fill_in_contribution_to_residuals(Vector<double>& residuals)
324 {
325 // Fill in the residuals of the Navier-Stokes equations
326 NavierStokesEquations<DIM>::fill_in_contribution_to_residuals(residuals);
327
328 // Fill in the residuals of the advection-diffusion eqautions
329 AdvectionDiffusionEquations<DIM>::fill_in_contribution_to_residuals(
330 residuals);
331 }
332
333
334//-----------Finite-difference the entire jacobian-----------------------
335//-----------------------------------------------------------------------
336#ifdef USE_FD_JACOBIAN_FOR_BUOYANT_Q_ELEMENT
337
338
339 /// Compute the element's residual vector and the Jacobian matrix.
340 /// Jacobian is computed by finite-differencing.
341 void fill_in_contribution_to_jacobian(Vector<double>& residuals,
342 DenseMatrix<double>& jacobian)
343 {
344 // This function computes the Jacobian by finite-differencing
345 FiniteElement::fill_in_contribution_to_jacobian(residuals, jacobian);
346 }
347
348 /// Add the element's contribution to its residuals vector,
349 /// jacobian matrix and mass matrix
351 Vector<double>& residuals,
352 DenseMatrix<double>& jacobian,
353 DenseMatrix<double>& mass_matrix)
354 {
355 // Call the standard (Broken) function
356 // which will prevent these elements from being used
357 // in eigenproblems until replaced.
358 FiniteElement::fill_in_contribution_to_jacobian_and_mass_matrix(
359 residuals, jacobian, mass_matrix);
360 }
361
362#else
363//--------Finite-difference off-diagonal-blocks in jacobian--------------
364//-----------------------------------------------------------------------
365#ifdef USE_OFF_DIAGONAL_FD_JACOBIAN_FOR_BUOYANT_Q_ELEMENT
366
367 /// Helper function to get the off-diagonal blocks of the Jacobian
368 /// matrix by finite differences
370 Vector<double>& residuals, DenseMatrix<double>& jacobian)
371 {
372 // Local storage for the index in the nodes at which the
373 // Navier-Stokes velocities are stored (we know that this should be 0,1,2)
374 unsigned u_nodal_nst[DIM];
375 for (unsigned i = 0; i < DIM; i++)
376 {
377 u_nodal_nst[i] = this->u_index_nst(i);
378 }
379
380 // Local storage for the index at which the temperature is stored
381 unsigned u_nodal_adv_diff = this->u_index_adv_diff();
382
383 // Find the total number of unknowns in the elements
384 unsigned n_dof = this->ndof();
385
386 // Temporary storage for residuals
387 Vector<double> newres(n_dof);
388
389 // Storage for local unknown and local equation
390 int local_unknown = 0, local_eqn = 0;
391
392 // Set the finite difference step
393 double fd_step = FiniteElement::Default_fd_jacobian_step;
394
395 // Find the number of nodes
396 unsigned n_node = this->nnode();
397
398 // Calculate the contribution of the Navier--Stokes velocities to the
399 // advection-diffusion equations
400
401 // Loop over the nodes
402 for (unsigned n = 0; n < n_node; n++)
403 {
404 // There are DIM values of the velocities
405 for (unsigned i = 0; i < DIM; i++)
406 {
407 // Get the local velocity equation number
408 local_unknown = this->nodal_local_eqn(n, u_nodal_nst[i]);
409
410 // If it's not pinned
411 if (local_unknown >= 0)
412 {
413 // Get a pointer to the velocity value
414 double* value_pt = this->node_pt(n)->value_pt(u_nodal_nst[i]);
415
416 // Save the old value
417 double old_var = *value_pt;
418
419 // Increment the value
420 *value_pt += fd_step;
421
422 // Get the altered advection-diffusion residuals.
423 // which must be done using fill_in_contribution because
424 // get_residuals will always return the full residuals
425 // because the appropriate fill_in function is overloaded above
426 for (unsigned m = 0; m < n_dof; m++)
427 {
428 newres[m] = 0.0;
429 }
430 AdvectionDiffusionEquations<DIM>::fill_in_contribution_to_residuals(
431 newres);
432
433 // AdvectionDiffusionEquations<DIM>::get_residuals(newres);
434
435 // Now fill in the Advection-Diffusion sub-block
436 // of the jacobian
437 for (unsigned m = 0; m < n_node; m++)
438 {
439 // Local equation for temperature
440 local_eqn = this->nodal_local_eqn(m, u_nodal_adv_diff);
441
442 // If it's not a boundary condition
443 if (local_eqn >= 0)
444 {
445 double sum =
446 (newres[local_eqn] - residuals[local_eqn]) / fd_step;
447 jacobian(local_eqn, local_unknown) = sum;
448 }
449 }
450
451 // Reset the nodal data
452 *value_pt = old_var;
453 }
454 }
455
456
457 // Calculate the contribution of the temperature to the Navier--Stokes
458 // equations
459 {
460 // Get the local equation number
461 local_unknown = this->nodal_local_eqn(n, u_nodal_adv_diff);
462
463 // If it's not pinned
464 if (local_unknown >= 0)
465 {
466 // Get a pointer to the concentration value
467 double* value_pt = this->node_pt(n)->value_pt(u_nodal_adv_diff);
468
469 // Save the old value
470 double old_var = *value_pt;
471
472 // Increment the value (Really need access)
473 *value_pt += fd_step;
474
475 // Get the altered Navier--Stokes residuals
476 // which must be done using fill_in_contribution because
477 // get_residuals will always return the full residuals
478 // because the appropriate fill_in function is overloaded above
479 for (unsigned m = 0; m < n_dof; m++)
480 {
481 newres[m] = 0.0;
482 }
483 NavierStokesEquations<DIM>::fill_in_contribution_to_residuals(
484 newres);
485
486 // NavierStokesEquations<DIM>::get_residuals(newres);
487
488 // Now fill in the Navier-Stokes sub-block
489 for (unsigned m = 0; m < n_node; m++)
490 {
491 // Loop over the fluid velocities
492 for (unsigned j = 0; j < DIM; j++)
493 {
494 // Local fluid equation
495 local_eqn = this->nodal_local_eqn(m, u_nodal_nst[j]);
496 if (local_eqn >= 0)
497 {
498 double sum =
499 (newres[local_eqn] - residuals[local_eqn]) / fd_step;
500 jacobian(local_eqn, local_unknown) = sum;
501 }
502 }
503 }
504
505 // Reset the nodal data
506 *value_pt = old_var;
507 }
508 }
509
510 } // End of loop over nodes
511 }
512
513
514 /// Compute the element's residual Vector and the Jacobian matrix.
515 /// Use finite-differencing only for the off-diagonal blocks.
516 void fill_in_contribution_to_jacobian(Vector<double>& residuals,
517 DenseMatrix<double>& jacobian)
518 {
519 // Calculate the Navier-Stokes contributions (diagonal block and
520 // residuals)
521 NavierStokesEquations<DIM>::fill_in_contribution_to_jacobian(residuals,
522 jacobian);
523
524 // Calculate the advection-diffusion contributions
525 //(diagonal block and residuals)
526 AdvectionDiffusionEquations<DIM>::fill_in_contribution_to_jacobian(
527 residuals, jacobian);
528
529 // We now fill in the off-diagonal (interaction) blocks by finite
530 // differences.
532 } // End of jacobian calculation
533
534
535 /// Add the element's contribution to its residuals vector,
536 /// jacobian matrix and mass matrix
538 Vector<double>& residuals,
539 DenseMatrix<double>& jacobian,
540 DenseMatrix<double>& mass_matrix)
541 {
542 // Get the analytic diagonal terms
543 NavierStokesEquations<
545 jacobian,
546 mass_matrix);
547
548 AdvectionDiffusionEquations<
550 jacobian,
551 mass_matrix);
552
553 // Now fill in the off-diagonal blocks
555 }
556
557
558 //--------------------Analytic jacobian---------------------------------
559//-----------------------------------------------------------------------
560#else
561
562 /// Helper function to get the off-diagonal blocks of the Jacobian
563 /// matrix analytically
565 Vector<double>& residuals, DenseMatrix<double>& jacobian)
566 {
567 // We now fill in the off-diagonal (interaction) blocks analytically
568 // This requires knowledge of exactly how the residuals are assembled
569 // within the parent elements and involves yet another loop over
570 // the integration points!
571
572 // Local storage for the index in the nodes at which the
573 // Navier-Stokes velocities are stored (we know that this should be 0,1,2)
574 unsigned u_nodal_nst[DIM];
575 for (unsigned i = 0; i < DIM; i++)
576 {
577 u_nodal_nst[i] = this->u_index_nst(i);
578 }
579
580 // Local storage for the index at which the temperature is stored
581 const unsigned u_nodal_adv_diff = this->u_index_adv_diff();
582
583 // Find out how many nodes there are
584 const unsigned n_node = this->nnode();
585
586 // Set up memory for the shape and test functions and their derivatives
587 Shape psif(n_node), testf(n_node);
588 DShape dpsifdx(n_node, DIM), dtestfdx(n_node, DIM);
589
590 // Number of integration points
591 const unsigned n_intpt = this->integral_pt()->nweight();
592
593 // Get Physical Variables from Element
594 double Ra = this->ra();
595 double Pe = this->pe();
596 Vector<double> gravity = this->g();
597
598 // Integers to store the local equations and unknowns
599 int local_eqn = 0, local_unknown = 0;
600
601 // Loop over the integration points
602 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
603 {
604 // Get the integral weight
605 double w = this->integral_pt()->weight(ipt);
606
607 // Call the derivatives of the shape and test functions
608 double J = this->dshape_and_dtest_eulerian_at_knot_nst(
609 ipt, psif, dpsifdx, testf, dtestfdx);
610
611 // Premultiply the weights and the Jacobian
612 double W = w * J;
613
614 // Calculate local values of temperature derivatives
615 // Allocate
616 Vector<double> interpolated_du_adv_diff_dx(DIM, 0.0);
617
618 // Loop over nodes
619 for (unsigned l = 0; l < n_node; l++)
620 {
621 // Get the nodal value
622 double u_value = this->raw_nodal_value(l, u_nodal_adv_diff);
623 // Loop over the derivative directions
624 for (unsigned j = 0; j < DIM; j++)
625 {
626 interpolated_du_adv_diff_dx[j] += u_value * dpsifdx(l, j);
627 }
628 }
629
630 // Assemble the jacobian terms
631
632 // Loop over the test functions
633 for (unsigned l = 0; l < n_node; l++)
634 {
635 // Assemble the contributions of the temperature to
636 // the Navier--Stokes equations (which arise through the buoyancy
637 // body-force term)
638
639 // Loop over the velocity components in the Navier--Stokes equtions
640 for (unsigned i = 0; i < DIM; i++)
641 {
642 // If it's not a boundary condition
643 local_eqn = this->nodal_local_eqn(l, u_nodal_nst[i]);
644 if (local_eqn >= 0)
645 {
646 // Loop over the velocity shape functions again
647 for (unsigned l2 = 0; l2 < n_node; l2++)
648 {
649 // We have only the temperature degree of freedom to consider
650 // If it's non-zero add in the contribution
651 local_unknown = this->nodal_local_eqn(l2, u_nodal_adv_diff);
652 if (local_unknown >= 0)
653 {
654 // Add contribution to jacobian matrix
655 jacobian(local_eqn, local_unknown) +=
656 -gravity[i] * psif(l2) * Ra * testf(l) * W;
657 }
658 }
659 }
660 }
661
662 // Assemble the contributions of the fluid velocities to the
663 // advection-diffusion equation for the temperature
664 {
665 local_eqn = this->nodal_local_eqn(l, u_nodal_adv_diff);
666 // If it's not pinned
667 if (local_eqn >= 0)
668 {
669 // Loop over the shape functions again
670 for (unsigned l2 = 0; l2 < n_node; l2++)
671 {
672 // Loop over the velocity degrees of freedom
673 for (unsigned i2 = 0; i2 < DIM; i2++)
674 {
675 // Get the local unknown
676 local_unknown = this->nodal_local_eqn(l2, u_nodal_nst[i2]);
677 // If it's not pinned
678 if (local_unknown >= 0)
679 {
680 // Add the contribution to the jacobian matrix
681 jacobian(local_eqn, local_unknown) -=
682 Pe * psif(l2) * interpolated_du_adv_diff_dx[i2] *
683 testf(l) * W;
684 }
685 }
686 }
687 }
688 }
689 }
690 }
691 }
692
693
694 /// Compute the element's residual Vector and the Jacobian matrix.
695 /// Use analytic expressions for the off-diagonal blocks
696 void fill_in_contribution_to_jacobian(Vector<double>& residuals,
697 DenseMatrix<double>& jacobian)
698 {
699 // Calculate the Navier-Stokes contributions (diagonal block and
700 // residuals)
701 NavierStokesEquations<DIM>::fill_in_contribution_to_jacobian(residuals,
702 jacobian);
703
704 // Calculate the advection-diffusion contributions
705 //(diagonal block and residuals)
706 AdvectionDiffusionEquations<DIM>::fill_in_contribution_to_jacobian(
707 residuals, jacobian);
708
709 // Fill in the off diagonal blocks analytically
711 }
712
713
714 /// Add the element's contribution to its residuals vector,
715 /// jacobian matrix and mass matrix
717 Vector<double>& residuals,
718 DenseMatrix<double>& jacobian,
719 DenseMatrix<double>& mass_matrix)
720 {
721 // Get the analytic diagonal terms for the jacobian and mass matrix
722 NavierStokesEquations<
724 jacobian,
725 mass_matrix);
726
727 AdvectionDiffusionEquations<
729 jacobian,
730 mass_matrix);
731
732 // Now fill in the off-diagonal blocks in the jacobian matrix
734 }
735
736#endif
737#endif
738 };
739
740 //=========================================================================
741 /// Set the default physical value to be zero
742 //=========================================================================
743 template<>
745 0.0;
746 template<>
748 0.0;
749
750
751 //=======================================================================
752 /// Face geometry of the 2D Buoyant Crouzeix_Raviart elements
753 //=======================================================================
754 template<unsigned int DIM>
755 class FaceGeometry<BuoyantQCrouzeixRaviartElement<DIM>>
756 : public virtual QElement<DIM - 1, 3>
757 {
758 public:
759 FaceGeometry() : QElement<DIM - 1, 3>() {}
760 };
761
762 //=======================================================================
763 /// Face geometry of the Face geometry of 2D Buoyant Crouzeix_Raviart elements
764 //=======================================================================
765 template<>
766 class FaceGeometry<FaceGeometry<BuoyantQCrouzeixRaviartElement<2>>>
767 : public virtual PointElement
768 {
769 public:
770 FaceGeometry() : PointElement() {}
771 };
772
773
774 /// //////////////////////////////////////////////////////////////////////
775 /// //////////////////////////////////////////////////////////////////////
776 /// //////////////////////////////////////////////////////////////////////
777
778
779 //============start_element_class============================================
780 /// A RefineableElement class that solves the
781 /// Boussinesq approximation of the Navier--Stokes
782 /// and energy equations by coupling two pre-existing classes.
783 /// The RefineableQAdvectionDiffusionElement
784 /// with bi-quadratic interpolation for the
785 /// scalar variable (temperature) and
786 /// RefineableQCrouzeixRaviartElement which solves the Navier--Stokes
787 /// equations using bi-quadratic interpolation for the velocities and a
788 /// discontinuous bi-linear interpolation for the pressure. Note that we are
789 /// free to choose the order in which we store the variables at the nodes. In
790 /// this case we choose to store the variables in the order fluid velocities
791 /// followed by temperature. We must, therefore, overload the function
792 /// AdvectionDiffusionEquations<DIM>::u_index_adv_diff() to indicate that
793 /// the temperature is stored at the DIM-th position not the 0-th. We do not
794 /// need to overload the corresponding function in the
795 /// NavierStokesEquations<DIM> class because the velocities are stored
796 /// first. Finally, we choose to use the flux-recovery calculation from the
797 /// fluid velocities to provide the error used in the mesh adaptation.
798 //==========================================================================
799 template<unsigned DIM>
801 : public virtual RefineableQAdvectionDiffusionElement<DIM, 3>,
802 public virtual RefineableQCrouzeixRaviartElement<DIM>
803 {
804 private:
805 /// Pointer to a new physical variable, the Rayleigh number
806 double* Ra_pt;
807
808 /// The static default value of the Rayleigh number
810
811 public:
812 /// Constructor: call the underlying constructors and
813 /// initialise the pointer to the Rayleigh number to address the default
814 /// value of 0.0
816 : RefineableQAdvectionDiffusionElement<DIM, 3>(),
817 RefineableQCrouzeixRaviartElement<DIM>()
818 {
820 }
821
822 /// The required number of values stored at the nodes is
823 /// the sum of the required values of the two single-physics elements. This
824 /// step is generic for any composed element of this type.
825 inline unsigned required_nvalue(const unsigned& n) const
826 {
827 return (RefineableQAdvectionDiffusionElement<DIM, 3>::required_nvalue(n) +
828 RefineableQCrouzeixRaviartElement<DIM>::required_nvalue(n));
829 }
830
831 /// Access function for the Rayleigh number (const version)
832 const double& ra() const
833 {
834 return *Ra_pt;
835 }
836
837 /// Access function for the pointer to the Rayleigh number
838 double*& ra_pt()
839 {
840 return Ra_pt;
841 }
842
843
844 /// Final override for disable ALE
846 {
847 // Disable ALE in both sets of equations
848 RefineableNavierStokesEquations<DIM>::disable_ALE();
849 RefineableAdvectionDiffusionEquations<DIM>::disable_ALE();
850 }
851
852 /// Final override for enable ALE
854 {
855 // Enable ALE in both sets of equations
856 RefineableNavierStokesEquations<DIM>::enable_ALE();
857 RefineableAdvectionDiffusionEquations<DIM>::enable_ALE();
858 }
859
860
861 /// Number of scalars/fields output by this element. Broken
862 /// virtual. Needs to be implemented for each new specific element type.
863 /// Temporary dummy
864 unsigned nscalar_paraview() const
865 {
866 throw OomphLibError(
867 "This function hasn't been implemented for this element",
868 OOMPH_CURRENT_FUNCTION,
869 OOMPH_EXCEPTION_LOCATION);
870
871 // Dummy unsigned
872 return 0;
873 }
874
875 /// Write values of the i-th scalar field at the plot points. Broken
876 /// virtual. Needs to be implemented for each new specific element type.
877 /// Temporary dummy
878 void scalar_value_paraview(std::ofstream& file_out,
879 const unsigned& i,
880 const unsigned& nplot) const
881 {
882 throw OomphLibError(
883 "This function hasn't been implemented for this element",
884 OOMPH_CURRENT_FUNCTION,
885 OOMPH_EXCEPTION_LOCATION);
886 }
887
888
889 /// Name of the i-th scalar field. Default implementation
890 /// returns V1 for the first one, V2 for the second etc. Can (should!) be
891 /// overloaded with more meaningful names.
892 std::string scalar_name_paraview(const unsigned& i) const
893 {
894 return "V" + StringConversion::to_string(i);
895 }
896
897 /// Overload the standard output function with the broken default
898 void output(std::ostream& outfile)
899 {
900 FiniteElement::output(outfile);
901 }
902
903 /// Output function:
904 /// x,y,u or x,y,z,u at Nplot^DIM plot points
905 void output(std::ostream& outfile, const unsigned& nplot)
906 {
907 // vector of local coordinates
908 Vector<double> s(DIM);
909
910 // Tecplot header info
911 outfile << this->tecplot_zone_string(nplot);
912
913 // Loop over plot points
914 unsigned num_plot_points = this->nplot_points(nplot);
915 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
916 {
917 // Get local coordinates of plot point
918 this->get_s_plot(iplot, nplot, s);
919
920 // Output the position of the plot point
921 for (unsigned i = 0; i < DIM; i++)
922 {
923 outfile << this->interpolated_x(s, i) << " ";
924 }
925
926 // Output the fluid velocities at the plot point
927 for (unsigned i = 0; i < DIM; i++)
928 {
929 outfile << this->interpolated_u_nst(s, i) << " ";
930 }
931
932 // Output the fluid pressure at the plot point
933 outfile << this->interpolated_p_nst(s) << " ";
934
935 // Output the temperature at the plot point
936 outfile << this->interpolated_u_adv_diff(s) << " " << std::endl;
937 }
938 outfile << std::endl;
939
940 // Write tecplot footer (e.g. FE connectivity lists)
941 this->write_tecplot_zone_footer(outfile, nplot);
942 }
943
944 /// C-style output function: Broken default
945 void output(FILE* file_pt)
946 {
947 FiniteElement::output(file_pt);
948 }
949
950 /// C-style output function: Broken default
951 void output(FILE* file_pt, const unsigned& n_plot)
952 {
953 FiniteElement::output(file_pt, n_plot);
954 }
955
956 /// Output function for an exact solution: Broken default
957 void output_fct(std::ostream& outfile,
958 const unsigned& Nplot,
959 FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
960 {
961 FiniteElement::output_fct(outfile, Nplot, exact_soln_pt);
962 }
963
964
965 /// Output function for a time-dependent exact solution.
966 /// Broken default
967 void output_fct(std::ostream& outfile,
968 const unsigned& Nplot,
969 const double& time,
970 FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
971 {
972 FiniteElement::output_fct(outfile, Nplot, time, exact_soln_pt);
973 }
974
975 /// Overload the index at which the temperature
976 /// variable is stored. We choose to store is after the fluid velocities.
977 unsigned u_index_adv_diff() const
978 {
979 return DIM;
980 }
981
982 /// Number of vertex nodes in the element is obtained from the
983 /// geometric element.
984 unsigned nvertex_node() const
985 {
986 return QElement<DIM, 3>::nvertex_node();
987 }
988
989 /// Pointer to the j-th vertex node in the element,
990 /// Call the geometric element's function.
991 Node* vertex_node_pt(const unsigned& j) const
992 {
993 return QElement<DIM, 3>::vertex_node_pt(j);
994 }
995
996 /// The total number of continously interpolated values is
997 /// DIM+1 (DIM fluid velocities and one temperature).
999 {
1000 return DIM + 1;
1001 }
1002
1003
1004 /// Get the continuously interpolated values at the local coordinate
1005 /// s. We choose to put the fluid velocities first, followed by the
1006 /// temperature.
1007 void get_interpolated_values(const Vector<double>& s,
1008 Vector<double>& values)
1009 {
1010 // Storage for the fluid velocities
1011 Vector<double> nst_values;
1012
1013 // Get the fluid velocities from the fluid element
1014 RefineableQCrouzeixRaviartElement<DIM>::get_interpolated_values(
1015 s, nst_values);
1016
1017 // Storage for the temperature
1018 Vector<double> advection_values;
1019
1020 // Get the temperature from the advection-diffusion element
1021 RefineableQAdvectionDiffusionElement<DIM, 3>::get_interpolated_values(
1022 s, advection_values);
1023
1024 // Add the fluid velocities to the values vector
1025 for (unsigned i = 0; i < DIM; i++)
1026 {
1027 values.push_back(nst_values[i]);
1028 }
1029
1030 // Add the concentration to the end
1031 values.push_back(advection_values[0]);
1032 }
1033
1034
1035 /// Get all continuously interpolated values at the local
1036 /// coordinate s at time level t (t=0: present; t>0: previous).
1037 /// We choose to put the fluid velocities first, followed by the
1038 /// temperature
1039 void get_interpolated_values(const unsigned& t,
1040 const Vector<double>& s,
1041 Vector<double>& values)
1042 {
1043 // Storage for the fluid velocities
1044 Vector<double> nst_values;
1045
1046 // Get the fluid velocities from the fluid element
1047 RefineableQCrouzeixRaviartElement<DIM>::get_interpolated_values(
1048 t, s, nst_values);
1049
1050 // Storage for the temperature
1051 Vector<double> advection_values;
1052
1053 // Get the temperature from the advection-diffusion element
1054 RefineableQAdvectionDiffusionElement<DIM, 3>::get_interpolated_values(
1055 s, advection_values);
1056
1057 // Add the fluid velocities to the values vector
1058 for (unsigned i = 0; i < DIM; i++)
1059 {
1060 values.push_back(nst_values[i]);
1061 }
1062
1063 // Add the concentration to the end
1064 values.push_back(advection_values[0]);
1065
1066 } // end of get_interpolated_values
1067
1068
1069 /// The additional hanging node information must be set up
1070 /// for both single-physics elements.
1072 {
1073 RefineableQCrouzeixRaviartElement<DIM>::further_setup_hanging_nodes();
1074 RefineableQAdvectionDiffusionElement<DIM,
1076 }
1077
1078
1079 /// Call the rebuild_from_sons functions for each of the
1080 /// constituent multi-physics elements.
1081 void rebuild_from_sons(Mesh*& mesh_pt)
1082 {
1083 RefineableQAdvectionDiffusionElement<DIM, 3>::rebuild_from_sons(mesh_pt);
1084 RefineableQCrouzeixRaviartElement<DIM>::rebuild_from_sons(mesh_pt);
1085 }
1086
1087
1088 /// Call the underlying single-physics element's further_build()
1089 /// functions and make sure that the pointer to the Rayleigh number
1090 /// is passed to the sons
1092 {
1093 RefineableQCrouzeixRaviartElement<DIM>::further_build();
1094 RefineableQAdvectionDiffusionElement<DIM, 3>::further_build();
1095
1096 // Cast the pointer to the father element to the specific
1097 // element type
1098 RefineableBuoyantQCrouzeixRaviartElement<DIM>* cast_father_element_pt =
1100 this->father_element_pt());
1101
1102 // Set the pointer to the Rayleigh number to be the same as that in
1103 // the father
1104 this->Ra_pt = cast_father_element_pt->ra_pt();
1105 } // end of further build
1106
1107
1108 /// The recovery order is that of the NavierStokes elements.
1110 {
1111 return RefineableQCrouzeixRaviartElement<DIM>::nrecovery_order();
1112 }
1113
1114 /// The number of Z2 flux terms is the same as that in
1115 /// the fluid element plus that in the advection-diffusion element
1117 {
1118 return (
1119 RefineableQCrouzeixRaviartElement<DIM>::num_Z2_flux_terms() +
1120 RefineableQAdvectionDiffusionElement<DIM, 3>::num_Z2_flux_terms());
1121 }
1122
1123
1124 /// Get the Z2 flux by concatenating the fluxes from the fluid and
1125 /// the advection diffusion elements.
1126 void get_Z2_flux(const Vector<double>& s, Vector<double>& flux)
1127 {
1128 // Find the number of fluid fluxes
1129 unsigned n_fluid_flux =
1130 RefineableQCrouzeixRaviartElement<DIM>::num_Z2_flux_terms();
1131
1132 // Fill in the first flux entries as the velocity entries
1133 RefineableQCrouzeixRaviartElement<DIM>::get_Z2_flux(s, flux);
1134
1135 // Find the number of temperature fluxes
1136 unsigned n_temp_flux =
1137 RefineableQAdvectionDiffusionElement<DIM, 3>::num_Z2_flux_terms();
1138 Vector<double> temp_flux(n_temp_flux);
1139
1140 // Get the temperature flux
1141 RefineableQAdvectionDiffusionElement<DIM, 3>::get_Z2_flux(s, temp_flux);
1142
1143 // Add the temperature flux to the end of the flux vector
1144 for (unsigned i = 0; i < n_temp_flux; i++)
1145 {
1146 flux[n_fluid_flux + i] = temp_flux[i];
1147 }
1148
1149 } // end of get_Z2_flux
1150
1151 /// The number of compound fluxes is two (one for the fluid and
1152 /// one for the temperature)
1154 {
1155 return 2;
1156 }
1157
1158 /// Fill in which flux components are associated with the fluid
1159 /// measure and which are associated with the temperature measure
1160 void get_Z2_compound_flux_indices(Vector<unsigned>& flux_index)
1161 {
1162 // Find the number of fluid fluxes
1163 unsigned n_fluid_flux =
1164 RefineableQCrouzeixRaviartElement<DIM>::num_Z2_flux_terms();
1165 // Find the number of temperature fluxes
1166 unsigned n_temp_flux =
1167 RefineableQAdvectionDiffusionElement<DIM, 3>::num_Z2_flux_terms();
1168
1169 // The fluid fluxes are first
1170 // The values of the flux_index vector are zero on entry, so we
1171 // could omit this line
1172 for (unsigned i = 0; i < n_fluid_flux; i++)
1173 {
1174 flux_index[i] = 0;
1175 }
1176
1177 // Set the temperature fluxes (the last set of fluxes
1178 for (unsigned i = 0; i < n_temp_flux; i++)
1179 {
1180 flux_index[n_fluid_flux + i] = 1;
1181 }
1182
1183 } // end of get_Z2_compound_flux_indices
1184
1185
1186 /// Validate against exact solution at given time
1187 /// Solution is provided via function pointer.
1188 /// Plot at a given number of plot points and compute L2 error
1189 /// and L2 norm of velocity solution over element
1190 /// Overload to broken default
1191 void compute_error(std::ostream& outfile,
1192 FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt,
1193 const double& time,
1194 double& error,
1195 double& norm)
1196 {
1197 FiniteElement::compute_error(outfile, exact_soln_pt, time, error, norm);
1198 }
1199
1200 /// Validate against exact solution.
1201 /// Solution is provided via function pointer.
1202 /// Plot at a given number of plot points and compute L2 error
1203 /// and L2 norm of velocity solution over element
1204 /// Overload to broken default.
1205 void compute_error(std::ostream& outfile,
1206 FiniteElement::SteadyExactSolutionFctPt exact_soln_pt,
1207 double& error,
1208 double& norm)
1209 {
1210 FiniteElement::compute_error(outfile, exact_soln_pt, error, norm);
1211 }
1212
1213 /// Overload the wind function in the advection-diffusion equations.
1214 /// This provides the coupling from the Navier--Stokes equations to the
1215 /// advection-diffusion equations because the wind is the fluid velocity.
1216 void get_wind_adv_diff(const unsigned& ipt,
1217 const Vector<double>& s,
1218 const Vector<double>& x,
1219 Vector<double>& wind) const
1220 {
1221 // The wind function is simply the velocity at the points
1222 this->interpolated_u_nst(s, wind);
1223 }
1224
1225
1226 /// Overload the body force in the navier-stokes equations
1227 /// This provides the coupling from the advection-diffusion equations
1228 /// to the Navier--Stokes equations, the body force is the
1229 /// temperature multiplied by the Rayleigh number acting in the
1230 /// direction opposite to gravity.
1231 void get_body_force_nst(const double& time,
1232 const unsigned& ipt,
1233 const Vector<double>& s,
1234 const Vector<double>& x,
1235 Vector<double>& result)
1236 {
1237 // Get the temperature
1238 const double interpolated_t = this->interpolated_u_adv_diff(s);
1239
1240 // Get vector that indicates the direction of gravity from
1241 // the Navier-Stokes equations
1242 Vector<double> gravity(NavierStokesEquations<DIM>::g());
1243
1244 // Temperature-dependent body force:
1245 for (unsigned i = 0; i < DIM; i++)
1246 {
1247 result[i] = -gravity[i] * interpolated_t * ra();
1248 }
1249 }
1250
1251 /// Fill in the constituent elements' contribution to the residual vector.
1252 void fill_in_contribution_to_residuals(Vector<double>& residuals)
1253 {
1254 // Call the residuals of the Navier-Stokes equations
1255 RefineableNavierStokesEquations<DIM>::fill_in_contribution_to_residuals(
1256 residuals);
1257
1258 // Call the residuals of the advection-diffusion equations
1259 RefineableAdvectionDiffusionEquations<
1261 }
1262
1263
1264 /// Compute the element's residual Vector and the jacobian matrix
1265 /// using full finite differences, the default implementation
1266 void fill_in_contribution_to_jacobian(Vector<double>& residuals,
1267 DenseMatrix<double>& jacobian)
1268 {
1269#ifdef USE_FD_JACOBIAN_FOR_REFINEABLE_BUOYANT_Q_ELEMENT
1270 FiniteElement::fill_in_contribution_to_jacobian(residuals, jacobian);
1271#else
1272 // Calculate the Navier-Stokes contributions (diagonal block and
1273 // residuals)
1274 RefineableNavierStokesEquations<DIM>::fill_in_contribution_to_jacobian(
1275 residuals, jacobian);
1276
1277 // Calculate the advection-diffusion contributions
1278 //(diagonal block and residuals)
1279 RefineableAdvectionDiffusionEquations<
1280 DIM>::fill_in_contribution_to_jacobian(residuals, jacobian);
1281
1282 // We now fill in the off-diagonal (interaction) blocks analytically
1283 this->fill_in_off_diagonal_jacobian_blocks_analytic(residuals, jacobian);
1284#endif
1285 } // End of jacobian calculation
1286
1287 /// Add the element's contribution to its residuals vector,
1288 /// jacobian matrix and mass matrix
1290 Vector<double>& residuals,
1291 DenseMatrix<double>& jacobian,
1292 DenseMatrix<double>& mass_matrix)
1293 {
1294 // Call the (broken) version in the base class
1295 FiniteElement::fill_in_contribution_to_jacobian_and_mass_matrix(
1296 residuals, jacobian, mass_matrix);
1297 }
1298
1299 /// Compute the contribution of the off-diagonal blocks
1300 /// analytically.
1302 Vector<double>& residuals, DenseMatrix<double>& jacobian)
1303 {
1304 // Perform another loop over the integration loops using the information
1305 // from the original elements' residual assembly loops to determine
1306 // the conributions to the jacobian
1307
1308 // Local storage for pointers to hang_info objects
1309 HangInfo *hang_info_pt = 0, *hang_info2_pt = 0;
1310
1311 // Local storage for the index in the nodes at which the
1312 // Navier-Stokes velocities are stored (we know that this should be 0,1,2)
1313 unsigned u_nodal_nst[DIM];
1314 for (unsigned i = 0; i < DIM; i++)
1315 {
1316 u_nodal_nst[i] = this->u_index_nst(i);
1317 }
1318
1319 // Local storage for the index at which the temperature is stored
1320 const unsigned u_nodal_adv_diff = this->u_index_adv_diff();
1321
1322 // Find out how many nodes there are
1323 const unsigned n_node = this->nnode();
1324
1325 // Set up memory for the shape and test functions and their derivatives
1326 Shape psif(n_node), testf(n_node);
1327 DShape dpsifdx(n_node, DIM), dtestfdx(n_node, DIM);
1328
1329 // Number of integration points
1330 const unsigned n_intpt = this->integral_pt()->nweight();
1331
1332 // Get Physical Variables from Element
1333 double Ra = this->ra();
1334 double Pe = this->pe();
1335 Vector<double> gravity = this->g();
1336
1337 // Integers to store the local equations and unknowns
1338 int local_eqn = 0, local_unknown = 0;
1339
1340 // Loop over the integration points
1341 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1342 {
1343 // Get the integral weight
1344 double w = this->integral_pt()->weight(ipt);
1345
1346 // Call the derivatives of the shape and test functions
1347 double J = this->dshape_and_dtest_eulerian_at_knot_nst(
1348 ipt, psif, dpsifdx, testf, dtestfdx);
1349
1350 // Premultiply the weights and the Jacobian
1351 double W = w * J;
1352
1353 // Calculate local values of temperature derivatives
1354 // Allocate
1355 Vector<double> interpolated_du_adv_diff_dx(DIM, 0.0);
1356
1357 // Loop over nodes
1358 for (unsigned l = 0; l < n_node; l++)
1359 {
1360 // Get the nodal value
1361 double u_value = this->nodal_value(l, u_nodal_adv_diff);
1362 // Loop over the derivative directions
1363 for (unsigned j = 0; j < DIM; j++)
1364 {
1365 interpolated_du_adv_diff_dx[j] += u_value * dpsifdx(l, j);
1366 }
1367 }
1368
1369 // Assemble the Jacobian terms
1370 //---------------------------
1371
1372 // Loop over the test functions/eqns
1373 for (unsigned l = 0; l < n_node; l++)
1374 {
1375 // Local variables to store the number of master nodes and
1376 // the weight associated with the shape function if the node is
1377 // hanging
1378 unsigned n_master = 1;
1379 double hang_weight = 1.0;
1380
1381 // Local bool (is the node hanging)
1382 bool is_node_hanging = this->node_pt(l)->is_hanging();
1383
1384 // If the node is hanging, get the number of master nodes
1385 if (is_node_hanging)
1386 {
1387 hang_info_pt = this->node_pt(l)->hanging_pt();
1388 n_master = hang_info_pt->nmaster();
1389 }
1390 // Otherwise there is just one master node, the node itself
1391 else
1392 {
1393 n_master = 1;
1394 }
1395
1396 // Loop over the master nodes
1397 for (unsigned m = 0; m < n_master; m++)
1398 {
1399 // If the node is hanging get weight from master node
1400 if (is_node_hanging)
1401 {
1402 // Get the hang weight from the master node
1403 hang_weight = hang_info_pt->master_weight(m);
1404 }
1405 else
1406 {
1407 // Node contributes with full weight
1408 hang_weight = 1.0;
1409 }
1410
1411
1412 // Assemble derivatives of Navier Stokes momentum w.r.t. temperature
1413 //-----------------------------------------------------------------
1414
1415 // Loop over velocity components for equations
1416 for (unsigned i = 0; i < DIM; i++)
1417 {
1418 // Get the equation number
1419 if (is_node_hanging)
1420 {
1421 // Get the equation number from the master node
1422 local_eqn = this->local_hang_eqn(
1423 hang_info_pt->master_node_pt(m), u_nodal_nst[i]);
1424 }
1425 else
1426 {
1427 // Local equation number
1428 local_eqn = this->nodal_local_eqn(l, u_nodal_nst[i]);
1429 }
1430
1431 if (local_eqn >= 0)
1432 {
1433 // Local variables to store the number of master nodes
1434 // and the weights associated with each hanging node
1435 unsigned n_master2 = 1;
1436 double hang_weight2 = 1.0;
1437
1438 // Loop over the nodes for the unknowns
1439 for (unsigned l2 = 0; l2 < n_node; l2++)
1440 {
1441 // Local bool (is the node hanging)
1442 bool is_node2_hanging = this->node_pt(l2)->is_hanging();
1443
1444 // If the node is hanging, get the number of master nodes
1445 if (is_node2_hanging)
1446 {
1447 hang_info2_pt = this->node_pt(l2)->hanging_pt();
1448 n_master2 = hang_info2_pt->nmaster();
1449 }
1450 // Otherwise there is one master node, the node itself
1451 else
1452 {
1453 n_master2 = 1;
1454 }
1455
1456 // Loop over the master nodes
1457 for (unsigned m2 = 0; m2 < n_master2; m2++)
1458 {
1459 if (is_node2_hanging)
1460 {
1461 // Read out the local unknown from the master node
1462 local_unknown = this->local_hang_eqn(
1463 hang_info2_pt->master_node_pt(m2), u_nodal_adv_diff);
1464 // Read out the hanging weight from the master node
1465 hang_weight2 = hang_info2_pt->master_weight(m2);
1466 }
1467 else
1468 {
1469 // The local unknown number comes from the node
1470 local_unknown =
1471 this->nodal_local_eqn(l2, u_nodal_adv_diff);
1472 // The hang weight is one
1473 hang_weight2 = 1.0;
1474 }
1475
1476 if (local_unknown >= 0)
1477 {
1478 // Add contribution to jacobian matrix
1479 jacobian(local_eqn, local_unknown) +=
1480 -gravity[i] * psif(l2) * Ra * testf(l) * W *
1481 hang_weight * hang_weight2;
1482 }
1483 }
1484 }
1485 }
1486 }
1487
1488
1489 // Assemble derivative of adv diff eqn w.r.t. fluid veloc
1490 //------------------------------------------------------
1491 {
1492 // Get the equation number
1493 if (is_node_hanging)
1494 {
1495 // Get the equation number from the master node
1496 local_eqn = this->local_hang_eqn(
1497 hang_info_pt->master_node_pt(m), u_nodal_adv_diff);
1498 }
1499 else
1500 {
1501 // Local equation number
1502 local_eqn = this->nodal_local_eqn(l, u_nodal_adv_diff);
1503 }
1504
1505 // If it's not pinned
1506 if (local_eqn >= 0)
1507 {
1508 // Local variables to store the number of master nodes
1509 // and the weights associated with each hanging node
1510 unsigned n_master2 = 1;
1511 double hang_weight2 = 1.0;
1512
1513 // Loop over the nodes for the unknowns
1514 for (unsigned l2 = 0; l2 < n_node; l2++)
1515 {
1516 // Local bool (is the node hanging)
1517 bool is_node2_hanging = this->node_pt(l2)->is_hanging();
1518
1519 // If the node is hanging, get the number of master nodes
1520 if (is_node2_hanging)
1521 {
1522 hang_info2_pt = this->node_pt(l2)->hanging_pt();
1523 n_master2 = hang_info2_pt->nmaster();
1524 }
1525 // Otherwise there is one master node, the node itself
1526 else
1527 {
1528 n_master2 = 1;
1529 }
1530
1531 // Loop over the master nodes
1532 for (unsigned m2 = 0; m2 < n_master2; m2++)
1533 {
1534 // If the node is hanging
1535 if (is_node2_hanging)
1536 {
1537 // Read out the hanging weight from the master node
1538 hang_weight2 = hang_info2_pt->master_weight(m2);
1539 }
1540 // If the node is not hanging
1541 else
1542 {
1543 // The hang weight is one
1544 hang_weight2 = 1.0;
1545 }
1546
1547 // Loop over the velocity degrees of freedom
1548 for (unsigned i2 = 0; i2 < DIM; i2++)
1549 {
1550 // If the node is hanging
1551 if (is_node2_hanging)
1552 {
1553 // Read out the local unknown from the master node
1554 local_unknown = this->local_hang_eqn(
1555 hang_info2_pt->master_node_pt(m2), u_nodal_nst[i2]);
1556 }
1557 else
1558 {
1559 // The local unknown number comes from the node
1560 local_unknown =
1561 this->nodal_local_eqn(l2, u_nodal_nst[i2]);
1562 }
1563
1564 // If it's not pinned
1565 if (local_unknown >= 0)
1566 {
1567 // Add the contribution to the jacobian matrix
1568 jacobian(local_eqn, local_unknown) -=
1569 Pe * psif(l2) * interpolated_du_adv_diff_dx[i2] *
1570 testf(l) * W * hang_weight * hang_weight2;
1571 }
1572 }
1573 }
1574 }
1575 }
1576 }
1577 }
1578 }
1579 }
1580 } // End of function
1581 };
1582
1583
1584 //===================================================================
1585 /// Set the default value of the Rayleigh number to be zero
1586 //===================================================================
1587 template<>
1588 double RefineableBuoyantQCrouzeixRaviartElement<
1590
1591 template<>
1594
1595
1596} // namespace oomph
1597
1598#endif
////////////////////////////////////////////////////////////////////// //////////////////////////////...
double * Ra_pt
Pointer to a private data member, the Rayleigh number.
double *& ra_pt()
Access function for the pointer to the Rayleigh number.
std::string scalar_name_paraview(const unsigned &i) const
Name of the i-th scalar field. Default implementation returns V1 for the first one,...
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Validate against exact solution. Solution is provided via function pointer. Plot at a given number of...
static double Default_Physical_Constant_Value
The static default value of the Rayleigh number.
void scalar_value_paraview(std::ofstream &file_out, const unsigned &i, const unsigned &nplot) const
Write values of the i-th scalar field at the plot points. Broken virtual. Needs to be implemented for...
void output(std::ostream &outfile)
Overload the standard output function with the broken default.
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function: Broken default.
unsigned nscalar_paraview() const
Number of scalars/fields output by this element. Broken virtual. Needs to be implemented for each new...
void disable_ALE()
Final override for disable ALE.
void compute_error(std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
Validate against exact solution at given time Solution is provided via function pointer....
void fill_in_off_diagonal_jacobian_blocks_by_fd(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Helper function to get the off-diagonal blocks of the Jacobian matrix by finite differences.
unsigned required_nvalue(const unsigned &n) const
The required number of values stored at the nodes is the sum of the required values of the two single...
void 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 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 output(std::ostream &outfile, const unsigned &nplot)
Output function: Output x, y, u, v, p, theta at Nplot^DIM plot points.
void unfix_pressure(const unsigned &p_dof)
Unpin p_dof-th pressure dof.
BuoyantQCrouzeixRaviartElement()
Constructor: call the underlying constructors and initialise the pointer to the Rayleigh number to po...
void output(FILE *file_pt)
C-style output function: Broken default.
void fill_in_off_diagonal_jacobian_blocks_analytic(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Helper function to get the off-diagonal blocks of the Jacobian matrix analytically.
const double & ra() const
Access function for the Rayleigh number (const version)
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Calculate the element's contribution to the residual vector. Recall that fill_in_* functions MUST NOT...
unsigned u_index_adv_diff() const
Overload the index at which the temperature variable is stored. We choose to store it after the fluid...
void output_fct(std::ostream &outfile, const unsigned &Nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output function for an exact solution: 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 output_fct(std::ostream &outfile, const unsigned &Nplot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
Output function for a time-dependent exact solution: Broken default.
void get_body_force_nst(const double &time, const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, Vector< double > &result)
Overload the body force in the Navier-Stokes equations This provides the coupling from the advection-...
void enable_ALE()
Final override for enable ALE.
////////////////////////////////////////////////////////////////////// //////////////////////////////...
void rebuild_from_sons(Mesh *&mesh_pt)
Call the rebuild_from_sons functions for each of the constituent multi-physics elements.
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 compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Validate against exact solution. Solution is provided via function pointer. Plot at a given number of...
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 using full finite differences,...
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Fill in the constituent elements' contribution to the residual vector.
RefineableBuoyantQCrouzeixRaviartElement()
Constructor: call the underlying constructors and initialise the pointer to the Rayleigh number to ad...
unsigned nrecovery_order()
The recovery order is that of the NavierStokes elements.
void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Get the Z2 flux by concatenating the fluxes from the fluid and the advection diffusion elements.
unsigned nscalar_paraview() const
Number of scalars/fields output by this element. Broken virtual. Needs to be implemented for each new...
static double Default_Physical_Constant_Value
The static default value of the Rayleigh number.
void output(std::ostream &outfile, const unsigned &nplot)
Output function: x,y,u or x,y,z,u at Nplot^DIM plot points.
unsigned ncompound_fluxes()
The number of compound fluxes is two (one for the fluid and one for the temperature)
const double & ra() const
Access function for the Rayleigh number (const version)
void further_build()
Call the underlying single-physics element's further_build() functions and make sure that the pointer...
unsigned num_Z2_flux_terms()
The number of Z2 flux terms is the same as that in the fluid element plus that in the advection-diffu...
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function: Broken default.
double * Ra_pt
Pointer to a new physical variable, the Rayleigh number.
unsigned required_nvalue(const unsigned &n) const
The required number of values stored at the nodes is the sum of the required values of the two single...
void get_body_force_nst(const double &time, const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, Vector< double > &result)
Overload the body force in the navier-stokes equations This provides the coupling from the advection-...
void fill_in_off_diagonal_jacobian_blocks_analytic(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Compute the contribution of the off-diagonal blocks analytically.
void enable_ALE()
Final override for enable ALE.
void disable_ALE()
Final override for disable ALE.
unsigned nvertex_node() const
Number of vertex nodes in the element is obtained from the geometric element.
void output(FILE *file_pt)
C-style output function: Broken default.
unsigned ncont_interpolated_values() const
The total number of continously interpolated values is DIM+1 (DIM fluid velocities and one temperatur...
void scalar_value_paraview(std::ofstream &file_out, const unsigned &i, const unsigned &nplot) const
Write values of the i-th scalar field at the plot points. Broken virtual. Needs to be implemented for...
void output_fct(std::ostream &outfile, const unsigned &Nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output function for an exact solution: Broken default.
void output_fct(std::ostream &outfile, const unsigned &Nplot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
Output function for a time-dependent exact solution. Broken default.
void get_interpolated_values(const Vector< double > &s, Vector< double > &values)
Get the continuously interpolated values at the local coordinate s. We choose to put the fluid veloci...
void further_setup_hanging_nodes()
The additional hanging node information must be set up for both single-physics elements.
double *& ra_pt()
Access function for the pointer to the Rayleigh number.
void output(std::ostream &outfile)
Overload the standard output function with the broken default.
unsigned u_index_adv_diff() const
Overload the index at which the temperature variable is stored. We choose to store is after the fluid...
void get_Z2_compound_flux_indices(Vector< unsigned > &flux_index)
Fill in which flux components are associated with the fluid measure and which are associated with the...
void compute_error(std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
Validate against exact solution at given time Solution is provided via function pointer....
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element, Call the geometric element's function.
std::string scalar_name_paraview(const unsigned &i) const
Name of the i-th scalar field. Default implementation returns V1 for the first one,...
void get_interpolated_values(const unsigned &t, const Vector< double > &s, Vector< double > &values)
Get all continuously interpolated values at the local coordinate s at time level t (t=0: present; t>0...
double Default_Physical_Constant_Value
Default value for physical constants.