dg_elements.cc
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// Non-inline member functions for discontinuous galerkin elements
27
28// oomph-lib includes
29#include "dg_elements.h"
30#include "shape.h"
31#include <iomanip>
32
33namespace oomph
34{
35 //===================================================================
36 /// Find pointers to neighbouring faces and the local coordinates
37 /// in those faces that correspond to the integration points in the
38 /// present face.
39 /// This is achieved by moving up to the bulk element and thence
40 /// the mesh which MUST have implemented a neighbour finding scheme
41 //==================================================================
43 const bool& add_neighbour_data_to_bulk)
44 {
45 // Cache the pointer to the bulk element
47 dynamic_cast<DGElement*>(this->bulk_element_pt());
48
49 // Find the number of points in the integration scheme
50 const unsigned n_intpt = integral_pt()->nweight();
51 // Resize the storage in the element
52 Neighbour_face_pt.resize(n_intpt);
53 Neighbour_local_coordinate.resize(n_intpt);
54
55 // If we are adding the neighbour data to the bulk element
56 // then resize this storage
57 if (add_neighbour_data_to_bulk)
58 {
59 Neighbour_external_data.resize(n_intpt);
60 }
61
62
63 // Get the dimension of the present element
64 const unsigned el_dim = this->dim();
65 // Local coordinate in the face element
66 Vector<double> s(el_dim);
67
68 // Get the dimension of the bulk element
69 const unsigned n_dim = bulk_element_pt->dim();
70 // Local coordinate in the bulk element
71 Vector<double> s_bulk(n_dim);
72
73 // Storage for the interpolation data in the neighbour
74 Vector<Data*> neighbour_data;
75
76 // Now loop over the integration points
77 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
78 {
79 // Get the local coordinate of the integration points
80 for (unsigned i = 0; i < el_dim; i++)
81 {
82 s[i] = integral_pt()->knot(ipt, i);
83 }
84
85 // Now get the bulk coordinate
86 this->get_local_coordinate_in_bulk(s, s_bulk);
87
88 // Now we find the neighbouring face via the bulk element's member
89 // function which calls the Mesh's member function
90 bulk_element_pt->get_neighbouring_face_and_local_coordinate(
91 this->face_index(),
92 s_bulk,
95
96 // If we are adding the external data to the bulk
97 if (add_neighbour_data_to_bulk)
98 {
99 // Throw the appropriate data from the neighbour face to the set
100 dynamic_cast<DGFaceElement*>(Neighbour_face_pt[ipt])
101 ->get_interpolation_data(neighbour_data);
102
103 // Find the number of data
104 unsigned n_neighbour_data = neighbour_data.size();
105 // Resize the storage accordingly
106 Neighbour_external_data.resize(n_neighbour_data);
107
108 // Add the data to the external data of the bulk element
109 for (unsigned n = 0; n < n_neighbour_data; n++)
110 {
112 bulk_element_pt->add_external_data(neighbour_data[n]);
113 }
114 }
115 }
116 }
117
118
119 //======================================================================
120 // Report the global coordinates corresponding to the integration points
121 // and the corresponding coordinates in the neighbouring faces
122 //======================================================================
124 {
125 // Find the number of nodes
126 const unsigned n_node = this->nnode();
127 // Allocate storage for the shape functions
128 Shape psi(n_node);
129
130 // Find the dimension of the problem
131 const unsigned dim = this->nodal_dimension();
132 // Storage for local coordinates in this face and its neighbour
133 Vector<double> x(dim), face_x(dim);
134
135 // Find the dimension of the element
136 const unsigned el_dim = this->dim();
137 // Storage for local coordinate
138 Vector<double> s(el_dim);
139
140 // Calculate the number of integration points from the array
141 const unsigned n_intpt = this->integral_pt()->nweight();
142 // Loop over the integration points
143 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
144 {
145 // Find the local coordinate at the knot point
146 for (unsigned i = 0; i < el_dim; i++)
147 {
148 s[i] = this->integral_pt()->knot(ipt, i);
149 }
150 // Get the value of the global coordinate in the present face
151 this->interpolated_x(s, x);
152
153 // Get the value of x in the neighbouring face
154 Neighbour_face_pt[ipt]->interpolated_x(Neighbour_local_coordinate[ipt],
155 face_x);
156
157 // Let's compare
158 oomph_info << "In Face In Neighbour\n";
159 for (unsigned i = 0; i < dim; i++)
160 {
161 if (i == 0)
162 {
163 oomph_info << "(";
164 }
165 else
166 {
167 oomph_info << ", ";
168 }
169 oomph_info << std::setw(5) << std::left << x[i];
170 }
171 oomph_info << ")";
172
173 oomph_info << " ";
174
175 for (unsigned i = 0; i < dim; i++)
176 {
177 if (i == 0)
178 {
179 oomph_info << "(";
180 }
181 else
182 {
183 oomph_info << ", ";
184 }
185 oomph_info << std::setw(5) << std::left << face_x[i];
186 }
187 oomph_info << ")";
188 oomph_info << std::endl;
189 }
190 }
191
192
193 //=====================================================================
194 /// Return the interpolated values of the unknown fluxes
195 //=====================================================================
197 {
198 // Find the number of nodes
199 const unsigned n_node = nnode();
200 // If there are no nodes then return immediately
201 if (n_node == 0)
202 {
203 return;
204 }
205
206 // Get the shape functions at the local coordinate
207 Shape psi(n_node);
208 this->shape(s, psi);
209
210 // Find the number of fluxes
211 const unsigned n_flux = this->required_nflux();
212
213 // Find the indices at which the local fluxes are stored
214 Vector<unsigned> flux_nodal_index(n_flux);
215 for (unsigned i = 0; i < n_flux; i++)
216 {
217 flux_nodal_index[i] = this->flux_index(i);
218 }
219 // Initialise the fluxes to zero
220 for (unsigned i = 0; i < n_flux; i++)
221 {
222 u[i] = 0.0;
223 }
224
225 // Loop over the nodes
226 for (unsigned n = 0; n < n_node; n++)
227 {
228 const double psi_ = psi[n];
229 for (unsigned i = 0; i < n_flux; i++)
230 {
231 u[i] += this->node_pt(n)->value(flux_nodal_index[i]) * psi_;
232 }
233 }
234 }
235
236 //=====================================================================
237 /// Add all the data that are used to interpolate the unknowns. This must
238 /// be consistent with the interpolated_u function above and the default
239 /// implementation assumes pure nodal interpolaton.
240 //=====================================================================
242 {
243 // Find the number of nodes
244 const unsigned n_node = nnode();
245 // Now resize the vector
246 interpolation_data.resize(n_node);
247
248 // If there are no nodes then return immediately
249 if (n_node == 0)
250 {
251 return;
252 }
253
254 // Otherwise loop over the nodes and add to the vector in order
255 for (unsigned n = 0; n < n_node; n++)
256 {
257 interpolation_data[n] = this->node_pt(n);
258 }
259 }
260
261 //====================================================================
262 /// Calculate the numerical flux at the knot point ipt. This is the
263 /// most general interface than can be overloaded if desired. The shape
264 /// functions at the knot point will be passed into this function.
265 //====================================================================
267 const Shape& psi,
268 Vector<double>& flux,
269 DenseMatrix<double>& dflux_du_int,
270 DenseMatrix<double>& dflux_du_ext,
271 unsigned flag)
272 {
273 // Find the number of nodes
274 const unsigned n_node = this->nnode();
275 // Find the nodal dimension
276 const unsigned nodal_dim = this->nodal_dimension();
277 // Number of fluxes
278 const unsigned n_flux = this->required_nflux();
279 // Find the indices at which the local fluxes are stored
280 Vector<unsigned> flux_nodal_index(n_flux);
281 for (unsigned i = 0; i < n_flux; i++)
282 {
283 flux_nodal_index[i] = this->flux_index(i);
284 }
285
286 // Calculate the local unknowns
287 Vector<double> interpolated_u(n_flux, 0.0);
288
289 // Loop over the shape functions
290 for (unsigned l = 0; l < n_node; l++)
291 {
292 // Cache the shape functions
293 const double psi_ = psi(l);
294 // Loop over the fluxes
295 for (unsigned i = 0; i < n_flux; i++)
296 {
297 // Calculate the velocity from the most recent nodal values
298 interpolated_u[i] += this->nodal_value(l, flux_nodal_index[i]) * psi_;
299 }
300 }
301
302 // Now calculate the outer unit normal Vector
303 Vector<double> interpolated_n(nodal_dim);
304 this->outer_unit_normal(ipt, interpolated_n);
305
306 // Get the pointer to the neighbour
307 DGFaceElement* neighbour_element_pt =
308 dynamic_cast<DGFaceElement*>(Neighbour_face_pt[ipt]);
309
310 // Get the neighbour's fluxes
311 Vector<double> interpolated_u_neigh(n_flux);
312
313 neighbour_element_pt->interpolated_u(Neighbour_local_coordinate[ipt],
314 interpolated_u_neigh);
315
316 // Call the "standard" numerical flux function
317 this->numerical_flux(
318 interpolated_n, interpolated_u, interpolated_u_neigh, flux);
319
320 // If we are calculating the jacobian add this term
321 if (flag && (flag < 3))
322 {
323 this->dnumerical_flux_du(interpolated_n,
325 interpolated_u_neigh,
326 dflux_du_int,
327 dflux_du_ext);
328 }
329 }
330
331 //===============================================================
332 /// Calculate the derivative of the
333 /// normal flux, which is the dot product of our
334 /// approximation to the flux with the outer unit normal,
335 /// with respect to the interior and exterior variables
336 /// Default is to use finite differences
337 //=============================================================
339 const Vector<double>& u_int,
340 const Vector<double>& u_ext,
341 DenseMatrix<double>& dflux_du_int,
342 DenseMatrix<double>& dflux_du_ext)
343 {
344 // Find the number of fluxes
345 const unsigned n_flux = this->required_nflux();
346
347 // Get a local copy of the unknowns
348 Vector<double> u_int_local = u_int;
349 Vector<double> u_ext_local = u_ext;
350
351 // Storage for incremented and decremented fluxes
352 Vector<double> flux_plus(n_flux), flux_minus(n_flux);
353
355
356 // Now loop over all the fluxes
357 for (unsigned n = 0; n < n_flux; n++)
358 {
359 // Increase internal value
360
361 // Store the old value
362 double old_var = u_int_local[n];
363 // Increment the value
364 u_int_local[n] += fd_step;
365 // Get the new values
366 this->numerical_flux(n_out, u_int_local, u_ext_local, flux_plus);
367
368 // Reset the value
369 u_int_local[n] = old_var;
370 // Decrement the value
371 u_int_local[n] -= fd_step;
372 // Get the new values
373 this->numerical_flux(n_out, u_int_local, u_ext_local, flux_minus);
374
375 // Assemble the column of the jacobian
376 for (unsigned m = 0; m < n_flux; m++)
377 {
378 dflux_du_int(m, n) = (flux_plus[m] - flux_minus[m]) / (2.0 * fd_step);
379 }
380
381 // Reset the value
382 u_int_local[n] = old_var;
383
384 // Increase external value
385
386 // Store the old value
387 old_var = u_ext_local[n];
388 // Increment the value
389 u_ext_local[n] += fd_step;
390 // Get the new values
391 this->numerical_flux(n_out, u_int_local, u_ext_local, flux_plus);
392
393 // Reset the value
394 u_ext_local[n] = old_var;
395 // Decrement the value
396 u_ext_local[n] -= fd_step;
397 // Get the new values
398 this->numerical_flux(n_out, u_int_local, u_ext_local, flux_minus);
399
400 // Assemble the column of the jacobian
401 for (unsigned m = 0; m < n_flux; m++)
402 {
403 dflux_du_ext(m, n) = (flux_plus[m] - flux_minus[m]) / (2.0 * fd_step);
404 }
405
406 // Reset the value
407 u_ext_local[n] = old_var;
408 }
409 }
410
411
412 //===================================================================
413 /// Calculate the integrated (numerical) flux out of the face and add
414 /// it to the residuals vector
415 //===================================================================
417 DenseMatrix<double>& jacobian,
418 unsigned flag)
419 {
420// Check that we have set up the coupling data if we are computing the jacobin
421#ifdef PARANOID
422 if (flag && (flag < 3))
423 {
424 if (Neighbour_external_data.size() == 0)
425 {
426 std::ostringstream error_stream;
427 error_stream
428 << "Coupling data between elements not included in jacobian\n"
429 << "You should call DGMesh::setup_face_neighbour_info(true) to "
430 "ensure\n"
431 << "that this information is included in the jacobian\n";
432
433 throw OomphLibError(
434 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
435 }
436 }
437#endif
438
439
440 // Find the number of nodes
441 const unsigned n_node = nnode();
442 // Dimension of the face
443 const unsigned el_dim = dim();
444 // Storage for the shape functions
445 Shape psi(n_node);
446
447 // Number of integration points
448 const unsigned n_intpt = this->integral_pt()->nweight();
449 // Number of fluxes
450 const unsigned n_flux = this->required_nflux();
451 // Find the indices at which the local fluxes are stored
452 Vector<unsigned> flux_nodal_index(n_flux);
453 for (unsigned i = 0; i < n_flux; i++)
454 {
455 flux_nodal_index[i] = this->flux_index(i);
456 }
457
458 // Cache the bulk element
459 DGElement* bulk_elem_pt = dynamic_cast<DGElement*>(this->bulk_element_pt());
460
461 // Storage for the flux and its derivatives
462 Vector<double> F(n_flux);
463 DenseMatrix<double> dF_du_int(n_flux, n_flux);
464 DenseMatrix<double> dF_du_ext(n_flux, n_flux);
465
466 // Loop over the integration points
467 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
468 {
469 // Get the integral weight
470 double W = this->integral_pt()->weight(ipt);
471 // Get the shape functions at the knot
472 this->shape_at_knot(ipt, psi);
473
474 // Calculate the Jacobian
475 // For a point element, it's one
476 double J = W;
477 // Otherwise calculate for the element
478 if (el_dim != 0)
479 {
480 J *= this->J_eulerian_at_knot(ipt);
481 }
482
483 // Now calculate the numerical flux (and derivatives)
484 this->numerical_flux_at_knot(ipt, psi, F, dF_du_int, dF_du_ext, flag);
485
486 // Limit if desired here
487
488 // Cache the pointer to the neighbour
489 DGFaceElement* neighbour_element_pt =
490 dynamic_cast<DGFaceElement*>(Neighbour_face_pt[ipt]);
491
492 // Finally we need to assemble the appropriate contributions
493 // to the residuals
494 for (unsigned l = 0; l < n_node; l++)
495 {
496 // Loop over the fluxes
497 for (unsigned i = 0; i < n_flux; i++)
498 {
499 // Get the local equation number in the bulk
500 int local_eqn = bulk_elem_pt->nodal_local_eqn(bulk_node_number(l),
501 flux_nodal_index[i]);
502
503 // If it's not a boundary condition
504 if (local_eqn >= 0)
505 {
506 // Add the flux multiplied by the shape function and the jacobian
507 residuals[local_eqn] -= psi(l) * F[i] * J;
508
509 // Add the jacobian contributions
510 if (flag)
511 {
512 // If we are assembling the jacobian
513 if (flag < 3)
514 {
515 // Loop over the internal nodes and fluxes again
516 for (unsigned l2 = 0; l2 < n_node; l2++)
517 {
518 for (unsigned i2 = 0; i2 < n_flux; i2++)
519 {
520 // Get the local unknown equation number in the bulk
521 int local_unknown = bulk_elem_pt->nodal_local_eqn(
522 bulk_node_number(l2), flux_nodal_index[i2]);
523
524 // If it's not a boundary condition
525 if (local_unknown >= 0)
526 {
527 // Add the flux multiplied by the shape function a
528 // nd the jacobian
529 jacobian(local_eqn, local_unknown) -=
530 psi(l) * dF_du_int(i, i2) * psi(l2) * J;
531 }
532 }
533 }
534
535 // How many nodes does it have
536 unsigned neigh_n_node = neighbour_element_pt->nnode();
537
538 // Get the flux indices from the neighbour
539 Vector<unsigned> neigh_flux_index(n_flux);
540 for (unsigned i2 = 0; i2 < n_flux; i2++)
541 {
542 neigh_flux_index[i2] = neighbour_element_pt->flux_index(i2);
543 }
544
545 // Loop over the neighbours nodes
546 for (unsigned l2 = 0; l2 < neigh_n_node; l2++)
547 {
548 // Loop over the fluxed
549 for (unsigned i2 = 0; i2 < n_flux; i2++)
550 {
551 // Get the local unknown equation number in the bulk
552 int local_unknown =
553 dynamic_cast<DGElement*>(this->bulk_element_pt())
555 flux_nodal_index[i2]);
556
557 // If it's not a boundary condition
558 if (local_unknown >= 0)
559 {
560 // Add the flux multiplied by the shape function
561 // and the jacobian
562 jacobian(local_eqn, local_unknown) -=
563 psi(l) * dF_du_ext(i, i2) * psi(l2) * J;
564 }
565 }
566 }
567 }
568 } // End of jacobian calculation
569 }
570 }
571 }
572 }
573 }
574
575 //========================================================================
576 /// Function that computes and stores the (inverse) mass matrix
577 //========================================================================
579 {
580 // Now let's assemble stuff
581 const unsigned n_dof = this->ndof();
582 // Allocate storage for the local mass matrix (if required)
583 if (M_pt == 0)
584 {
586 }
587
588 // Resize and initialise the vector that will holds the residuals
589 Vector<double> dummy(n_dof, 0.0);
590
591 // Resize the mass matrix
592 M_pt->resize(n_dof, n_dof);
593 // Initialise the entries to zero
594 M_pt->initialise(0.0);
595 // Get the local mass matrix and residuals
597
598 // Now invert the mass matrix it will always be small
599 // This can possibly be streamlined (for example in spectral
600 // elements the mass matrix is diagonal)
601 M_pt->ludecompose();
602
603 // The mass matrix has been computed
605 }
606
607
608 //============================================================================
609 /// Function that returns the current value of the residuals
610 /// multiplied by the inverse mass matrix (virtual so that it can be
611 /// overloaded specific elements in which time and memory saving tricks can be
612 /// applied)
613 //============================================================================
615 Vector<double>& minv_res)
616 {
617 // If there are external data this is not going to work
618 if (nexternal_data() > 0)
619 {
620 std::ostringstream error_stream;
621 error_stream
622 << "Cannot use a discontinuous formulation for the mass matrix when\n"
623 << "there are external data.\n "
624 << "Do not call Problem::enable_discontinuous_formulation()\n";
625
626 throw OomphLibError(
627 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
628 }
629
630 // Now let's assemble stuff
631 const unsigned n_dof = this->ndof();
632 // Allocate storage for the local mass matrix (if required)
633 if (M_pt == 0)
634 {
636 }
637
638 // Resize and initialise the vector that will holds the residuals
639 minv_res.resize(n_dof);
640 for (unsigned n = 0; n < n_dof; n++)
641 {
642 minv_res[n] = 0.0;
643 }
644
645 // If we are recycling the mass matrix
647 {
648 // Get the residuals
649 this->fill_in_contribution_to_residuals(minv_res);
650 }
651 // Otherwise
652 else
653 {
654 // Resize the mass matrix
655 M_pt->resize(n_dof, n_dof);
656 // Initialise the entries to zero
657 M_pt->initialise(0.0);
658 // Get the local mass matrix and residuals
660
661 // Now invert the mass matrix it will always be small
662 // This can possibly be streamlined (for example in spectral
663 // elements the mass matrix is diagonal)
664 M_pt->ludecompose();
665
666 // The mass matrix has been computed
668 }
669
670 // Always do the backsubstitution
671 M_pt->lubksub(minv_res);
672 }
673
674
676 const int& face_index,
677 const Vector<double>& s,
678 FaceElement*& face_element_pt,
679 Vector<double>& s_face)
680 {
681 DG_mesh_pt->neighbour_finder(this, face_index, s, face_element_pt, s_face);
682 }
683
684
685 /// Limit the slope within an element
686 void DGElement::slope_limit(SlopeLimiter* const& slope_limiter_pt)
687 {
688 // Firstly find the dimension
689 const unsigned n_dim = this->dim();
690 // Find the number of fluxes
691 const unsigned n_flux = this->required_nflux();
692
693 switch (n_dim)
694 {
695 // One dimensional (easy-ish) case
696 case 1:
697 {
698 // Storage for the element and its neighbours
699 Vector<DGElement*> required_element_pt(3);
700 required_element_pt[0] = this;
701
702 // Get the pointer to the element on the left
703 required_element_pt[1] = dynamic_cast<DGElement*>(
704 dynamic_cast<DGFaceElement*>(this->face_element_pt(0))
705 ->neighbour_face_pt(0)
706 ->bulk_element_pt());
707 // Get the pointer to the element on the right
708 required_element_pt[2] = dynamic_cast<DGElement*>(
709 dynamic_cast<DGFaceElement*>(this->face_element_pt(1))
710 ->neighbour_face_pt(0)
711 ->bulk_element_pt());
712
713 // Loop over the fluxed
714 for (unsigned i = 0; i < n_flux; i++)
715 {
716 // Call our limiter, which will take as it's arguments, the current
717 // element and the required neighbours
718 slope_limiter_pt->limit(i, required_element_pt);
719 }
720 }
721 break;
722
723 default:
724 {
725 std::ostringstream error_stream;
726 error_stream << "Slope limiting is not implemented for this dimension: "
727 << n_dim;
728 throw OomphLibError(
729 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
730 }
731 }
732 }
733
734 double DGMesh::FaceTolerance = 1.0e-10;
735
736
737 //====================================================
738 /// Helper minmod function
739 //====================================================
741 {
742 const unsigned n_arg = args.size();
743 // If no arguments return zero
744 if (n_arg == 0)
745 {
746 return 0.0;
747 }
748
749 // Initialise the sign from the sign of the first entry
750 int sign = 0;
751 if (args[0] < 0.0)
752 {
753 sign = -1;
754 }
755 else if (args[0] > 0.0)
756 {
757 sign = 1;
758 }
759 else
760 {
761 return 0.0;
762 }
763
764 // Initialise the minimum value
765 double min = std::fabs(args[0]);
766
767 // Now loop over the rest of the values
768 for (unsigned i = 1; i < n_arg; i++)
769 {
770 if (sign == 1)
771 {
772 if (args[i] < 0.0)
773 {
774 return 0.0;
775 }
776 else if (args[i] < min)
777 {
778 min = args[i];
779 }
780 }
781 else if (sign == -1)
782 {
783 if (args[i] > 0.0)
784 {
785 return 0.0;
786 }
787 else if (std::fabs(args[i]) < min)
788 {
789 min = std::fabs(args[i]);
790 }
791 }
792 }
793
794 // Make sure to return the sign multiplied by the minimum value
795 return sign * min;
796 }
797
798
799 /// Modified minmod limiter to fix behaviour in smooth regions
800 double MinModLimiter::minmodB(Vector<double>& args, const double& h)
801 {
802 const unsigned n_arg = args.size();
803 // If no arguments return zero
804 if (n_arg == 0)
805 {
806 return 0.0;
807 }
808
809 // Modification to fix extrema
810 if (std::fabs(args[0]) < this->M * h * h)
811 {
812 return args[0];
813 }
814 // Otherwise just return the usual minmod
815 else
816 {
817 return minmod(args);
818 }
819 }
820
821 /// Implement the limiter function for the basic MinModlimiter
822 void MinModLimiter::limit(const unsigned& i,
823 const Vector<DGElement*>& required_element_pt)
824 {
825 // Set the tolerance
826 const double tol = 1.0e-16;
827
828 // Find the geometric parameters of the element
829 const unsigned n_node = required_element_pt[0]->nnode();
830 const double x_l = required_element_pt[0]->node_pt(0)->x(0);
831 const double x_r = required_element_pt[0]->node_pt(n_node - 1)->x(0);
832 const double h = x_r - x_l;
833 const double x0 = 0.5 * (x_l + x_r);
834 // Find the average value
835 const double u_av = required_element_pt[0]->average_value(i);
836
837 // Storage for the gradients to the minmod function
838 Vector<double> arg;
839 arg.reserve(3);
840
841 // Add the approximation for the element's left flux
842 arg.push_back(u_av - required_element_pt[0]->node_pt(0)->value(i));
843
844 // If there is a left element add the approximate gradient
845 // to the argument vector
846 if (required_element_pt[1] != required_element_pt[0])
847 {
848 arg.push_back(u_av - required_element_pt[1]->average_value(i));
849 }
850 // If there is a right element add the approximate gradient
851 // to the argument vector
852 if (required_element_pt[2] != required_element_pt[0])
853 {
854 arg.push_back(required_element_pt[2]->average_value(i) - u_av);
855 }
856
857 // Set the left value
858 const double u_l = u_av - this->minmod(arg);
859
860 // Now replace the first term in the argument list with the
861 // approximation for the element's right flux
862 arg.front() = required_element_pt[0]->node_pt(n_node - 1)->value(i) - u_av;
863
864 // Set the right value
865 const double u_r = u_av + this->minmod(arg);
866
867 // If the basic limited values are different from
868 // the unlimited values then limit
869 if ((std::fabs(u_l - required_element_pt[0]->node_pt(0)->value(i)) > tol) &&
870 (std::fabs(
871 u_r - required_element_pt[0]->node_pt(n_node - 1)->value(i)) > tol))
872 {
873 // Find the centre of the element on the left
874 const double x0_l =
875 0.5 * (required_element_pt[1]
876 ->node_pt(required_element_pt[1]->nnode() - 1)
877 ->x(0) +
878 required_element_pt[1]->node_pt(0)->x(0));
879 // Find the centre of the element on the right
880 const double x0_r =
881 0.5 * (required_element_pt[2]
882 ->node_pt(required_element_pt[2]->nnode() - 1)
883 ->x(0) +
884 required_element_pt[2]->node_pt(0)->x(0));
885
886 // Clear the argument list and reserve its size to
887 arg.clear();
888 arg.reserve(3);
889
890 // Approximate the gradient over the whole cell
891 arg.push_back((required_element_pt[0]->node_pt(n_node - 1)->value(i) -
892 required_element_pt[0]->node_pt(0)->value(i)) /
893 h);
894
895 // Adjust the estimates for the gradient calculated from neighbouring
896 // averages for the straight min-mod and MUSCL limiters
897 double gradient_factor = 0.5;
898 if (MUSCL)
899 {
900 gradient_factor = 1.0;
901 }
902
903 // If there is a left element, form the gradient
904 if (required_element_pt[0] != required_element_pt[1])
905 {
906 arg.push_back((u_av - required_element_pt[1]->average_value(i)) /
907 (gradient_factor * (x0 - x0_l)));
908 }
909
910 // If there is a right element, form the gradient
911 if (required_element_pt[0] != required_element_pt[2])
912 {
913 arg.push_back((required_element_pt[2]->average_value(i) - u_av) /
914 (gradient_factor * (x0_r - x0)));
915 }
916
917 // Calculate the limited gradient of these three
918 double limited_gradient = this->minmodB(arg, h);
919
920 // Loop over the nodes and limit
921 for (unsigned n = 0; n < n_node; n++)
922 {
923 double x = required_element_pt[0]->node_pt(n)->x(0) - x0;
924 required_element_pt[0]->node_pt(n)->set_value(
925 i, u_av + x * limited_gradient);
926 }
927 }
928 }
929
930} // namespace oomph
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
A Base class for DGElements.
Definition: dg_elements.h:161
DGMesh * DG_mesh_pt
Pointer to Mesh, which will be responsible for the neighbour finding.
Definition: dg_elements.h:171
void pre_compute_mass_matrix()
Function that computes and stores the (inverse) mass matrix.
Definition: dg_elements.cc:578
void slope_limit(SlopeLimiter *const &slope_limiter_pt)
Limit the slope within the element.
Definition: dg_elements.cc:686
virtual unsigned required_nflux()
Set the number of flux components.
Definition: dg_elements.h:193
bool Mass_matrix_has_been_computed
Boolean flag to indicate whether the mass matrix has been computed.
Definition: dg_elements.h:186
DGFaceElement * face_element_pt(const unsigned &i)
Access function for the faces.
Definition: dg_elements.h:366
bool Mass_matrix_reuse_is_enabled
Boolean flag to indicate whether to reuse the mass matrix.
Definition: dg_elements.h:182
void get_neighbouring_face_and_local_coordinate(const int &face_index, const Vector< double > &s, FaceElement *&face_element_pt, Vector< double > &s_face)
Return the neighbour info.
Definition: dg_elements.cc:675
virtual void get_inverse_mass_matrix_times_residuals(Vector< double > &minv_res)
Function that returns the current value of the residuals multiplied by the inverse mass matrix (virtu...
Definition: dg_elements.cc:614
DenseDoubleMatrix * M_pt
Pointer to storage for a mass matrix that can be recycled if desired.
Definition: dg_elements.h:175
Base class for Discontinuous Galerkin Faces. These are responsible for calculating the normal fluxes ...
Definition: dg_elements.h:50
virtual void numerical_flux(const Vector< double > &n_out, const Vector< double > &u_int, const Vector< double > &u_ext, Vector< double > &flux)
Calculate the normal flux, which is the dot product of our approximation to the flux with the outer u...
Definition: dg_elements.h:121
virtual void get_interpolation_data(Vector< Data * > &interpolation_data)
Get the data that are used to interpolate the unkowns in the element. These must be returned in order...
Definition: dg_elements.cc:241
void setup_neighbour_info(const bool &add_neighbour_data_to_bulk)
Setup information from the neighbouring face, return a set of nodes (as data) in the neighour....
Definition: dg_elements.cc:42
void report_info()
Output information about the present element and its neighbour.
Definition: dg_elements.cc:123
Vector< FaceElement * > Neighbour_face_pt
Vector of neighbouring face elements at the integration points.
Definition: dg_elements.h:52
void add_flux_contributions(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
Add the contribution from integrating the numerical flux.
Definition: dg_elements.cc:416
Vector< Vector< unsigned > > Neighbour_external_data
Vector of the vectors that will store the number of the bulk external data that correspond to the dof...
Definition: dg_elements.h:63
virtual void dnumerical_flux_du(const Vector< double > &n_out, const Vector< double > &u_int, const Vector< double > &u_ext, DenseMatrix< double > &dflux_du_int, DenseMatrix< double > &dflux_du_ext)
Calculate the derivative of the normal flux, which is the dot product of our approximation to the flu...
Definition: dg_elements.cc:338
virtual void interpolated_u(const Vector< double > &s, Vector< double > &f)
Return the interpolated values of the unknown fluxes.
Definition: dg_elements.cc:196
virtual void numerical_flux_at_knot(const unsigned &ipt, const Shape &psi, Vector< double > &flux, DenseMatrix< double > &dflux_du_int, DenseMatrix< double > &dflux_du_ext, unsigned flag)
Calculate the normal numerical flux at the integration point. This is the most general interface that...
Definition: dg_elements.cc:266
Vector< Vector< double > > Neighbour_local_coordinate
Vector of neighbouring local coordinates at the integration points.
Definition: dg_elements.h:55
virtual unsigned flux_index(const unsigned &i) const
Return the index at which the i-th unknown flux is stored.
Definition: dg_elements.h:68
virtual unsigned required_nflux()
Set the number of flux components.
Definition: dg_elements.h:74
virtual void neighbour_finder(FiniteElement *const &bulk_element_pt, const int &face_index, const Vector< double > &s_bulk, FaceElement *&face_element_pt, Vector< double > &s_face)
Definition: dg_elements.h:470
static double FaceTolerance
Definition: dg_elements.h:464
Class of matrices containing doubles, and stored as a DenseMatrix<double>, but with solving functiona...
Definition: matrices.h:1271
virtual void lubksub(DoubleVector &rhs)
LU backsubstitution.
Definition: matrices.cc:202
virtual void ludecompose()
LU decomposition using DenseLU (default linea solver)
Definition: matrices.cc:192
void initialise(const T &val)
Initialize all values in the matrix to val.
Definition: matrices.h:514
void resize(const unsigned long &n)
Resize to a square nxn matrix; any values already present will be transfered.
Definition: matrices.h:498
FaceElements are elements that coincide with the faces of higher-dimensional "bulk" elements....
Definition: elements.h:4338
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
Definition: elements.h:4626
void outer_unit_normal(const Vector< double > &s, Vector< double > &unit_normal) const
Compute outer unit normal at the specified local coordinate.
Definition: elements.cc:6006
double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s. Overloaded to get information from bulk...
Definition: elements.h:4528
unsigned & bulk_node_number(const unsigned &n)
Return the bulk node number that corresponds to the n-th local node number.
Definition: elements.h:4825
FiniteElement *& bulk_element_pt()
Pointer to higher-dimensional "bulk" element.
Definition: elements.h:4735
double J_eulerian_at_knot(const unsigned &ipt) const
Return the Jacobian of the mapping from local to global coordinates at the ipt-th integration point O...
Definition: elements.cc:5328
void get_local_coordinate_in_bulk(const Vector< double > &s, Vector< double > &s_bulk) const
Calculate the vector of local coordinate in the bulk element given the local coordinates in this Face...
Definition: elements.cc:6384
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 shape(const Vector< double > &s, Shape &psi) const =0
Calculate the geometric shape functions at local coordinate s. This function must be overloaded for e...
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Return the local equation number corresponding to the i-th value at the n-th local node.
Definition: elements.h: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
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
unsigned nodal_dimension() const
Return the required Eulerian dimension of the nodes in this element.
Definition: elements.h:2484
virtual void shape_at_knot(const unsigned &ipt, Shape &psi) const
Return the geometric shape function at the ipt-th integration point.
Definition: elements.cc:3220
unsigned nexternal_data() const
Return the number of external data objects.
Definition: elements.h:829
virtual void fill_in_contribution_to_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &mass_matrix)
Add the elemental contribution to the mass matrix matrix. and the residuals vector....
Definition: elements.cc:1292
static double Default_fd_jacobian_step
Double used for the default finite difference step in elemental jacobian calculations.
Definition: elements.h:1198
virtual void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the elemental contribution to the residuals vector. Note that this function will NOT initialise t...
Definition: elements.h:357
unsigned ndof() const
Return the number of equations/dofs in the element.
Definition: elements.h:835
int external_local_eqn(const unsigned &i, const unsigned &j)
Return the local equation number corresponding to the j-th value stored at the i-th external data.
Definition: elements.h:311
unsigned add_external_data(Data *const &data_pt, const bool &fd=true)
Add a (pointer to an) external data object to the element and return its index (i....
Definition: elements.cc:307
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
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.
double minmod(Vector< double > &args)
The basic minmod function.
Definition: dg_elements.cc:740
bool MUSCL
Boolean flag to indicate a MUSCL or straight min-mod limiter.
Definition: dg_elements.h:550
double minmodB(Vector< double > &args, const double &h)
The modified minmod function.
Definition: dg_elements.cc:800
void limit(const unsigned &i, const Vector< DGElement * > &required_element_pt)
The limit function.
Definition: dg_elements.cc:822
double M
Constant that is used in the modified min-mod function to provide better properties at extrema.
Definition: dg_elements.h:547
double value(const unsigned &i) const
Return i-th value (dofs or pinned) at this node either directly or via hanging node representation....
Definition: nodes.cc:2408
An OomphLibError object which should be thrown when an run-time error is encountered....
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
Base class for slope limiters.
Definition: dg_elements.h:525
virtual void limit(const unsigned &i, const Vector< DGElement * > &required_element_pt)
Basic function.
Definition: dg_elements.h:534
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...