periodic_orbit_handler.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 
31 #include "periodic_orbit_handler.h"
32 
33 namespace oomph
34 {
35 
37  std::ostream &outfile,
38  const unsigned &n_plot)
39  {
40  //Find out how many nodes there are
41  const unsigned n_node = nnode();
42 
43  //Set up memory for the shape and test functions
44  Shape psi(n_node);
45  DShape dpsidt(n_node,1);
46 
47  Vector<double> s(1);
48 
49  //Cast the element
50  PeriodicOrbitBaseElement* const base_el_pt =
51  dynamic_cast<PeriodicOrbitBaseElement*>(elem_pt);
52 
53  const double inverse_timescale = this->omega();
54 
55  //Loop over the plot point
56  for(unsigned i=0;i<n_plot;i++)
57  {
58  //Get the coordinate
59  s[0] = -1.0 + 2.0*i/((double)(n_plot-1));
60 
61  (void)dshape_eulerian(s,psi,dpsidt);
62 
63  //Sort out the timestepper weights and the time in here
64  //Calculate the time
65  double interpolated_time = 0.0;
66  for(unsigned l=0;l<n_node;l++)
67  {
68  interpolated_time += raw_nodal_position(l,0)*psi(l);
69  }
70 
71 
72  //Need to multiply by period and the set global time
73  this->time_pt()->time() = interpolated_time;
74 
75  //Set the weights of the timestepper
76  this->set_timestepper_weights(psi,dpsidt);
77 
78  base_el_pt->spacetime_output(outfile,n_plot,
79  interpolated_time/inverse_timescale);
80  }
81  }
82 
84  PeriodicOrbitAssemblyHandlerBase* const &assembly_handler_pt,
85  GeneralisedElement* const &elem_pt,
86  Vector<double> &residuals, DenseMatrix<double> &jacobian,
87  const unsigned& flag)
88  {
89  //A simple integration loop
90  //Find out how many nodes there are
91  const unsigned n_node = nnode();
92 
93  //Set up memory for the shape and test functions
94  Shape psi(n_node), test(n_node);
95  DShape dpsidt(n_node,1), dtestdt(n_node,1);
96 
97  //Set the value of n_intpt
98  const unsigned n_intpt = integral_pt()->nweight();
99 
100  //Integers to store the local equation and unknown numbers
101  int local_eqn=0, local_unknown=0;
102 
103  //Storage for the underlying element's residuals
104  const unsigned n_elem_dof = elem_pt->ndof();
105  Vector<double> el_residuals(n_elem_dof);
106  DenseMatrix<double> el_mass;
107  DenseMatrix<double> el_jacobian;
108 
109  if(flag)
110  {
111  el_mass.resize(n_elem_dof,n_elem_dof);
112  el_jacobian.resize(n_elem_dof,n_elem_dof);
113  }
114 
115  //Storage for the current value of the unkowns
116  Vector<double> u(n_elem_dof);
117  //and the derivatives
118  Vector<double> du_dt(n_elem_dof);
119  //And the previous derivative
120  Vector<double> du_dt_old(n_elem_dof);
121  //Storage for the inner product matrix
122  DenseMatrix<double> inner_product(n_elem_dof);
123  //Cast the element
124  PeriodicOrbitBaseElement* const base_el_pt =
125  dynamic_cast<PeriodicOrbitBaseElement*>(elem_pt);
126  //Get the inner product matrix
127  base_el_pt->get_inner_product_matrix(inner_product);
128 
129  //Storage for "all" unknowns
130  Vector<double> all_current_unknowns;
131  Vector<double> all_previous_unknowns;
132 
133  //Get the value of omega
134  const double inverse_timescale = this->omega();
135 
136  //Loop over the integration points
137  for(unsigned ipt=0;ipt<n_intpt;ipt++)
138  {
139  //Get the integral weight
140  double w = integral_pt()->weight(ipt);
141 
142  //Call the derivatives of the shape and test functions
143  double J = dshape_and_dtest_eulerian_at_knot_orbit(ipt,psi,dpsidt,
144  test,dtestdt);
145 
146  //Premultiply the weights and the Jacobian
147  double W = w*J;
148 
149  //Sort out the timestepper weights and the time in here
150  //Calculate the time
151  double interpolated_time = 0.0;
152  for(unsigned l=0;l<n_node;l++)
153  {
154  interpolated_time += raw_nodal_position(l,0)*psi(l);
155  }
156 
157 
158  //Need to multiply by period and the set global time
159  this->time_pt()->time() = interpolated_time/inverse_timescale;
160  //Set the weights of the timestepper
161  this->set_timestepper_weights(psi,dpsidt);
162 
163  //Get the residuals of the element or residuals and jacobian
164  if(flag)
165  {
166  elem_pt->get_jacobian_and_mass_matrix(el_residuals,el_jacobian,
167  el_mass);
168  }
169  else
170  {
171  elem_pt->get_residuals(el_residuals);
172  }
173 
174  //Get the current value of the unknown time derivatives
175  base_el_pt->get_non_external_ddofs_dt(du_dt);
176 
177  //Multiply the elemental residuals by the appropriate test functions
178  for(unsigned l=0;l<n_node;l++)
179  {
180  //Get the local equation number
181  local_eqn = nodal_local_eqn(l,0);
182  //If it's not a boundary condition (it will never be)
183  if(local_eqn >= 0)
184  {
185  //Work out the offset which is the GLOBAL time unknown
186  //multiplied by the number of elemental unknowns
187  unsigned offset = this->eqn_number(local_eqn)*n_elem_dof;
188  //Add to the appropriate residuals
189  for(unsigned i=0;i<n_elem_dof;i++)
190  {
191  residuals[i + offset] += el_residuals[i]*psi(l)*W;
192  }
193 
194 
195  //Add the jacobian contributions
196  if(flag)
197  {
198  //The form of the equations is -M du/dt + R = 0
199 
200  //Now add the contribution to the jacobian
201  for(unsigned l2=0;l2<n_node;l2++)
202  {
203  local_unknown = nodal_local_eqn(l2,0);
204  if(local_unknown >= 0)
205  {
206  //Work out the second offset
207  unsigned offset2 = this->eqn_number(local_unknown)*n_elem_dof;
208  //Add to the appropriate jacobian terms
209  for(unsigned i=0;i<n_elem_dof;i++)
210  {
211  for(unsigned j=0;j<n_elem_dof;j++)
212  {
213  //Add in the Jacobian terms
214  jacobian(i + offset, j+ offset2) +=
215  el_jacobian(i,j)*psi(l2)*psi(l)*W;
216 
217  //Add the time derivative terms,
218  jacobian(i + offset, j + offset2) -= el_mass(i,j)*
219  dpsidt(l2,0)*inverse_timescale*psi(l)*W;
220  }
221  }
222  }
223  }
224 
225 
226  //Add the variation of the period
227  for(unsigned i=0;i<n_elem_dof;i++)
228  {
229  for(unsigned j=0;j<n_elem_dof;j++)
230  {
231  jacobian(i + offset, Ntstorage*n_elem_dof)
232  -= el_mass(i,j)*(du_dt[j]/inverse_timescale)*psi(l)*W;
233  }
234  }
235  } //End of Jacobian flag
236  }
237  }
238 
239  //Sort out the phase condition
240 
241  //Get the current value of the unknowns
242  base_el_pt->get_non_external_dofs(u);
243 
244  //Now get the unknowns required by the assembly handler
245  //i.e. including all values throughout the period for backup
246  assembly_handler_pt->get_dofs_for_element(elem_pt,
247  all_current_unknowns);
248  //Get the previous values as stored by the assembly handler
249  assembly_handler_pt->get_previous_dofs_for_element(elem_pt,
250  all_previous_unknowns);
251 
252  //Now set the elemental values to the previous values
253  assembly_handler_pt->set_dofs_for_element(elem_pt,all_previous_unknowns);
254  //Get the previous time derivatives
255  base_el_pt->get_non_external_ddofs_dt(du_dt_old);
256  //Reset the element's values to the current
257  assembly_handler_pt->set_dofs_for_element(elem_pt,all_current_unknowns);
258 
259 
260  //Assemble the inner product
261  double sum = 0.0;
262  for(unsigned i=0;i<n_elem_dof;i++)
263  {
264  for(unsigned j=0;j<n_elem_dof;j++)
265  {
266  sum += u[i]*inner_product(i,j)*du_dt_old[j];
267  }
268  }
269 
270  //Add to the residuals
271  residuals[Ntstorage*n_elem_dof] += sum*W;
272 
273  //Sort out the jacobian
274  if(flag)
275  {
276  //Loop over the unknown time points
277  for(unsigned l2=0;l2<n_node;l2++)
278  {
279  //Get the local unknown
280  local_unknown = nodal_local_eqn(l2,0);
281  if(local_unknown >= 0)
282  {
283  //Work out the offset
284  unsigned offset2 = this->eqn_number(local_unknown)*n_elem_dof;
285  //Now add in the appropriate jacobian terms
286  for(unsigned i2=0;i2<n_elem_dof;i2++)
287  {
288  double sum2 = 0.0;
289  for(unsigned j=0;j<n_elem_dof;j++)
290  {
291  sum2 += inner_product(i2,j)*du_dt_old[j];
292  }
293  jacobian(Ntstorage*n_elem_dof,i2 + offset2) += psi(l2)*sum2*W;
294  }
295  }
296  }
297  } //End of jacobian calculation
298 
299  } //End of loop over time period integration points
300  }
301 
302 
303 
304 }
virtual void get_inner_product_matrix(DenseMatrix< double > &inner_product)
Get the inner product matrix.
A Generalised Element class.
Definition: elements.h:76
virtual void get_residuals(Vector< double > &residuals)
Calculate the vector of residuals of the equations in the element. By default initialise the vector t...
Definition: elements.h:975
double omega()
Return the frequency.
double dshape_eulerian(const Vector< double > &s, Shape &psi, DShape &dpsidx) const
Compute the geometric shape functions and also first derivatives w.r.t. global coordinates at local c...
Definition: elements.cc:3225
cstr elem_len * i
Definition: cfortran.h:607
virtual void get_non_external_dofs(Vector< double > &u)
Interface to get the current value of all (internal and shared) unknowns.
unsigned long eqn_number(const unsigned &ieqn_local) const
Return the global equation number corresponding to the ieqn_local-th local equation number...
Definition: elements.h:709
virtual double dshape_and_dtest_eulerian_at_knot_orbit(const unsigned &ipt, Shape &psi, DShape &dpsidt, Shape &test, DShape &dtestdt) const =0
Shape/test functions and derivs w.r.t. to global coords at integration point ipt; return Jacobian of ...
double raw_nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. Do not use the hanging node representation. NOTE: Moved to cc file because of a possible compiler bug in gcc (yes, really!). The move to the cc file avoids inlining which appears to cause problems (only) when compiled with gcc and -O3; offensive "illegal read" is in optimised-out section of code and data that is allegedly illegal is readily readable (by other means) just before this function is called so I can't really see how we could possibly be responsible for this...
Definition: elements.cc:1657
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
Time *& time_pt()
Retun the pointer to the global time.
void orbit_output(GeneralisedElement *const &elem_pt, std::ostream &outfile, const unsigned &n_plot)
virtual void get_non_external_ddofs_dt(Vector< double > &du_dt)
static char t char * s
Definition: cfortran.h:572
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1900
virtual void get_dofs_for_element(GeneralisedElement *const elem_pt, Vector< double > &dofs)=0
virtual void get_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
Calculate the residuals and jacobian and elemental "mass" matrix, the matrix that multiplies the time...
Definition: elements.h:1010
double & time()
Return the current value of the continuous time.
Definition: timesteppers.h:130
virtual void set_dofs_for_element(GeneralisedElement *const elem_pt, Vector< double > const &dofs)=0
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
virtual void spacetime_output(std::ostream &outilfe, const unsigned &Nplot, const double &time=0.0)
unsigned ndof() const
Return the number of equations/dofs in the element.
Definition: elements.h:834
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
void set_timestepper_weights(const Shape &psi, const DShape &dpsidt)
Set the timestepper weights.
void fill_in_generic_residual_contribution_orbit(PeriodicOrbitAssemblyHandlerBase *const &assembly_handler_pt, GeneralisedElement *const &elem_pt, Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
The routine that actually does all the work!
virtual void get_previous_dofs_for_element(GeneralisedElement *const elem_pt, Vector< double > &dofs)=0
void resize(const unsigned long &n)
Definition: matrices.h:511