flux_transport_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 function of the flux transport elements class
31 
33 #include "../generic/shape.h"
34 #include "../generic/timesteppers.h"
35 
36 namespace oomph
37 {
38  //======================================================================
39  ///Return the derivatives of the flux components as functions of the
40  ///unknowns. The default implementation is to use second-order
41  ///finite-differences, but these can be overloaded in specific equations
42  //========================================================================
43 template<unsigned DIM>
46 {
47  //Find the number of fluxes
48  const unsigned n_flux = nflux();
49  //Local copy of the unknowns
50  Vector<double> u_local = u;
51  //Finite differences
52  DenseMatrix<double> F(n_flux,DIM), F_plus(n_flux,DIM), F_minus(n_flux,DIM);
53  const double fd_step = GeneralisedElement::Default_fd_jacobian_step;
54  //Now loop over all the fluxes
55  for(unsigned p=0;p<n_flux;p++)
56  {
57  //Store the old value
58  const double old_var = u_local[p];
59  //Increment the value
60  u_local[p] += fd_step;
61  //Get the new values
62  flux(u_local,F_plus);
63 
64  //Reset the value
65  u_local[p] = old_var;
66  //Decrement the value
67  u_local[p] -= fd_step;
68  //Get the new values
69  flux(u_local,F_minus);
70 
71  //Assemble the entries of the jacobian
72  for(unsigned r=0;r<n_flux;r++)
73  {
74  for(unsigned j=0;j<DIM;j++)
75  {
76  df_du(r,j,p) = (F_plus(r,j) - F_minus(r,j))/(2.0*fd_step);
77  }
78  }
79 
80  //Reset the value
81  u_local[p] = old_var;
82  }
83 }
84 
85  //=====================================================================
86  ///\short Compute the residuals for the flux transport equations;
87  /// flag=1 compute jacobian as well
88  /// flag=2 compute mass matrix and jacobian as well
89  /// flag=3 compute mass matrix as well
90  //=======================================================================
91  template<unsigned DIM>
94  Vector<double> &residuals, DenseMatrix<double> &jacobian,
95  DenseMatrix<double> &mass_matrix, unsigned flag)
96  {
97  //Find the number of fluxes
98  const unsigned n_flux = this->nflux();
99  //Find the number of nodes
100  const unsigned n_node = this->nnode();
101  //Storage for the shape function and derivatives of shape function
102  Shape psi(n_node), test(n_node);
103  DShape dpsidx(n_node,DIM), dtestdx(n_node,DIM);
104 
105  //Cache the nodal indices at which the unknowns are stored
106  unsigned u_nodal_index[n_flux];
107  for(unsigned i=0;i<n_flux;i++)
108  {u_nodal_index[i] = this->u_index_flux_transport(i);}
109 
110  //Find the number of integration points
111  unsigned n_intpt = this->integral_pt()->nweight();
112 
113  //Integers to store the local equations and unknowns
114  int local_eqn=0, local_unknown = 0;
115 
116  //Loop over the integration points
117  for(unsigned ipt=0;ipt<n_intpt;ipt++)
118  {
119  //Get the shape functions at the knot
120  double J =
121  this->dshape_and_dtest_eulerian_at_knot_flux_transport(ipt,psi,dpsidx,
122  test,dtestdx);
123 
124  //Get the integral weight
125  double W = this->integral_pt()->weight(ipt)*J;
126 
127  //Storage for the local unknowns
128  Vector<double> interpolated_u(n_flux,0.0);
129  Vector<double> interpolated_dudt(n_flux,0.0);
130 
131  //Loop over the shape functions
132  for(unsigned l=0;l<n_node;l++)
133  {
134  //Locally cache the shape function
135  const double psi_ = psi(l);
136  for(unsigned i=0;i<n_flux;i++)
137  {
138  //Calculate the velocity and tangent vector
139  interpolated_u[i] += this->nodal_value(l,u_nodal_index[i])*psi_;
140  interpolated_dudt[i] += this->du_dt_flux_transport(l,i)*psi_;
141  }
142  }
143 
144  //Calculate the flux at the integration point
145  DenseMatrix<double> F(n_flux,DIM);
146  this->flux(interpolated_u,F);
147 
148  RankThreeTensor<double> dF_du(n_flux,DIM,n_flux);
149  if((flag) && (flag!=3))
150  {this->dflux_du(interpolated_u,dF_du);}
151 
152  //We need to assemble the contributions to the volume integral
153  for(unsigned l=0;l<n_node;l++)
154  {
155  //Loop over the unknowns
156  for(unsigned i=0;i<n_flux;i++)
157  {
158  //Get the local equation number
159  local_eqn = this->nodal_local_eqn(l,i);
160 
161  //if it's not a boundary condition
162  if(local_eqn >= 0)
163  {
164  //Add the time derivatives
165  residuals[local_eqn] -= interpolated_dudt[i]*test(l)*W;
166 
167 
168  //Calculate the flux dot product
169  double flux_dot = 0.0;
170  for(unsigned k=0;k<DIM;k++)
171  {flux_dot += F(i,k)*dtestdx(l,k);}
172 
173  //Add the contribution to the residuals
174  residuals[local_eqn] += flux_dot*W;
175 
176  //Now worry about the jacobian and mass matrix terms
177  if(flag)
178  {
179  //If we are assembling the jacobian
180  if(flag < 3)
181  {
182  //Loop over the shape functions again
183  for(unsigned l2=0;l2<n_node;l2++)
184  {
185  //Loop over the unknowns again
186  for(unsigned i2=0;i2<n_flux;i2++)
187  {
188  //Get the local unknowns
189  local_unknown = this->nodal_local_eqn(l2,i2);
190  //If not pinned
191  if(local_unknown >= 0)
192  {
193  //Add the time derivative terms
194  if(i2==i)
195  {
196  jacobian(local_eqn,local_unknown) -=
197  node_pt(l2)->time_stepper_pt()->weight(1,0)*
198  psi(l2)*test(l)*W;
199 
200  //Add the mass matrix if we are computing
201  //it and the jacobian
202  if(flag==2)
203  {
204  mass_matrix(local_eqn,local_unknown) +=
205  psi(l2)*test(l)*W;
206  }
207  }
208 
209  //Add the flux derivative terms
210  double dflux_dot = 0.0;
211  for(unsigned k=0;k<DIM;k++)
212  {
213  dflux_dot += dF_du(i,k,i2)*dtestdx(l,k);
214  }
215  jacobian(local_eqn,local_unknown) +=
216  dflux_dot*psi(l2)*W;
217  }
218  }
219  }
220  }
221  //End of jacobian assembly, here we are just assembling the
222  //mass matrix
223  else
224  {
225  for(unsigned l2=0;l2<n_node;l2++)
226  {
227  local_unknown = this->nodal_local_eqn(l2,i);
228  if(local_unknown >= 0)
229  {
230  mass_matrix(local_eqn,local_unknown) += psi(l2)*test(l)*W;
231  }
232  }
233  }
234  }
235  }
236  }
237  }
238  }
239  }
240 
241  //==================================================================
242  ///Return the i-th unknown at the local coordinate s
243  //==================================================================
244 template<unsigned DIM>
246  const Vector<double> &s, const unsigned &i)
247  {
248  //Find the number of nodes
249  const unsigned n_node = this->nnode();
250  //Get the shape functions at the local coordinate
251  Shape psi(n_node);
252  this->shape(s,psi);
253 
254  //Now interpolate each unknown
255  double u = 0.0;
256 
257  //Loop over the nodes
258  for(unsigned n=0;n<n_node;n++)
259  {
260  const double psi_ = psi[n];
261  u += this->nodal_value(n,u_index_flux_transport(i))*psi_;
262  }
263  return u;
264  }
265 
266  //======================================================================
267  /// \short i-th component of du/dt at local node n.
268  /// Uses suitably interpolated value for hanging nodes.
269  //======================================================================
270  template<unsigned DIM>
272  du_dt_flux_transport( const unsigned &n, const unsigned &i) const
273  {
274  // Get the data's timestepper
275  TimeStepper* time_stepper_pt = this->node_pt(n)->time_stepper_pt();
276 
277  //Initialise dudt
278  double dudt=0.0;
279 
280  //Loop over the timesteps, if there is a non Steady timestepper
281  if (!time_stepper_pt->is_steady())
282  {
283  //Find the index at which the dof is stored
284  const unsigned u_nodal_index = this->u_index_flux_transport(i);
285 
286  // Number of timsteps (past & present)
287  const unsigned n_time = time_stepper_pt->ntstorage();
288  // Loop over the timesteps
289  for(unsigned t=0;t<n_time;t++)
290  {
291  dudt+=time_stepper_pt->weight(1,t)*nodal_value(t,n,u_nodal_index);
292 
293  }
294  }
295 
296  return dudt;
297  }
298 
299  //==================================================================
300  ///Calculate the averages of the flux variables
301  //=================================================================
302  template<unsigned DIM>
304  double* &average_value)
305  {
306  //Find the number of fluxes
307  const unsigned n_flux = this->nflux();
308  //Resize the memory if necessary
309  if(average_value==0) {average_value = new double[n_flux+1];}
310 
311  //Initialise the averages to zero
312  for(unsigned i=0;i<n_flux+1;i++) {average_value[i] = 0.0;}
313 
314  //Find the number of nodes
315  const unsigned n_node = this->nnode();
316  //Storage for the shape functions
317  Shape psi(n_node);
318  DShape dpsidx(n_node,DIM);
319 
320  //Cache the nodal indices at which the unknowns are stored
321  unsigned u_nodal_index[n_flux];
322  for(unsigned i=0;i<n_flux;i++)
323  {u_nodal_index[i] = this->u_index_flux_transport(i);}
324 
325  //Find the number of integration points
326  unsigned n_intpt = this->integral_pt()->nweight();
327 
328  //Loop over the integration points
329  for(unsigned ipt=0;ipt<n_intpt;ipt++)
330  {
331  //Get the shape functions and derivatives at the knot
332  double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
333 
334  //Get the integral weight
335  double W = this->integral_pt()->weight(ipt)*J;
336 
337  //Storage for the local unknowns
338  Vector<double> interpolated_u(n_flux,0.0);
339 
340  //Loop over the shape functions
341  for(unsigned l=0;l<n_node;l++)
342  {
343  //Locally cache the shape function
344  const double psi_ = psi(l);
345  for(unsigned i=0;i<n_flux;i++)
346  {
347  //Calculate the velocity and tangent vector
348  interpolated_u[i] += this->nodal_value(l,u_nodal_index[i])*psi_;
349  }
350  }
351 
352  average_value[n_flux] += W;
353  //Loop over the values
354  for(unsigned i=0;i<n_flux;i++)
355  {
356  average_value[i] += interpolated_u[i]*W;
357  }
358  }
359 
360  //Divide the results by the size of the element
361  for(unsigned i=0;i<n_flux;i++)
362  {
363  average_value[i] /= average_value[n_flux];
364  }
365  }
366 
367 
368 
369 
370  //====================================================================
371  ///Output function, print the values of all unknowns
372  //==================================================================
373  template<unsigned DIM>
375  output(std::ostream &outfile, const unsigned &nplot)
376  {
377  //Find the number of fluxes
378  const unsigned n_flux = this->nflux();
379 
380  //Vector of local coordinates
381  Vector<double> s(DIM);
382 
383  // Tecplot header info
384  outfile << tecplot_zone_string(nplot);
385 
386  // Loop over plot points
387  unsigned num_plot_points=nplot_points(nplot);
388  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
389  {
390 
391  // Get local coordinates of plot point
392  get_s_plot(iplot,nplot,s);
393 
394  // Coordinates
395  for(unsigned i=0;i<DIM;i++)
396  {
397  outfile << interpolated_x(s,i) << " ";
398  }
399 
400  // Values
401  for(unsigned i=0;i<n_flux;i++)
402  {
403  outfile << interpolated_u_flux_transport(s,i) << " ";
404  }
405 
406  outfile << std::endl;
407  }
408  outfile << std::endl;
409 
410  // Write tecplot footer (e.g. FE connectivity lists)
411  write_tecplot_zone_footer(outfile,nplot);
412 
413 }
414 
415  template class FluxTransportEquations<1>;
416  template class FluxTransportEquations<2>;
417  template class FluxTransportEquations<3>;
418 
419 }
void output(std::ostream &outfile)
Output the element data — typically the values at the nodes in a format suitable for post-processing...
cstr elem_len * i
Definition: cfortran.h:607
virtual double weight(const unsigned &i, const unsigned &j) const
Access function for j-th weight for the i-th derivative.
Definition: timesteppers.h:557
char t
Definition: cfortran.h:572
bool is_steady() const
Flag to indicate if a timestepper has been made steady (possibly temporarily to switch off time-depen...
Definition: timesteppers.h:375
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:562
A Rank 3 Tensor class.
Definition: matrices.h:1337
static char t char * s
Definition: cfortran.h:572
void shape(const double &s, double *Psi)
Definition for 1D Lagrange shape functions. The value of all the shape functions at the local coordin...
Definition: shape.h:549
double du_dt_flux_transport(const unsigned &n, const unsigned &i) const
i-th component of du/dt at local node n. Uses suitably interpolated value for hanging nodes...
double interpolated_u_flux_transport(const Vector< double > &s, const unsigned &i)
Return the i-th unknown at the local coordinate s.
virtual void fill_in_generic_residual_contribution_flux_transport(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
Compute the residuals for the Navier–Stokes equations; flag=1(or 0): do (or don't) compute the Jacobi...
void calculate_element_averages(double *&average_values)
Compute the average values of the fluxes.
virtual void dflux_du(const Vector< double > &u, RankThreeTensor< double > &df_du)
Return the flux derivatives as a funciton of the unknowns Default finite-different implementation...
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition: timesteppers.h:219
static double Default_fd_jacobian_step
Double used for the default finite difference step in elemental jacobian calculations.
Definition: elements.h:1164