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