ode_elements.h
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// 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 #ifndef OOMPH_ODE_ELEMENTS_H
31 #define OOMPH_ODE_ELEMENTS_H
32 
33 #include "../generic/oomph_definitions.h"
34 #include "../generic/oomph_utilities.h"
35 
36 #include "../generic/matrices.h"
37 #include "../generic/Vector.h"
38 #include "../generic/elements.h"
39 #include "../generic/timesteppers.h"
40 
41 namespace oomph
42 {
43 
44  /// Element for integrating an initial value ODE
46  {
47 
48  public:
49 
50  /// Default constructor: null any pointers
52  {
53  Exact_solution_pt = 0;
54 
55  Use_fd_jacobian = false;
56  }
57 
58  /// Constructor: Pass time stepper and a solution function pointer, then
59  /// build the element.
60  ODEElement(TimeStepper* time_stepper_pt,
61  SolutionFunctorBase* exact_solution_pt)
62  {
63  build(time_stepper_pt, exact_solution_pt);
64  }
65 
66  /// Store pointers, create internal data.
67  void build(TimeStepper* time_stepper_pt,
68  SolutionFunctorBase* exact_solution_pt)
69  {
70  Exact_solution_pt = exact_solution_pt;
71 
72  Vector<double> exact = this->exact_solution(0);
73  unsigned nvalue = exact.size();
74 
75  add_internal_data(new Data(time_stepper_pt, nvalue));
76 
77  Use_fd_jacobian = false;
78  }
79 
80 
81  virtual ~ODEElement() {}
82 
83  unsigned nvalue() const
84  {
85  return internal_data_pt(0)->nvalue();
86  }
87 
88  /// Get residuals
90  {
91  // Get pointer to one-and-only internal data object
92  Data* dat_pt = internal_data_pt(0);
93 
94  // Get it's values
95  Vector<double> u(nvalue(), 0.0);
96  dat_pt->value(u);
97 
98  // Get time stepper
99  TimeStepper* time_stepper_pt = dat_pt->time_stepper_pt();
100 
101  // Get continuous time
102  double t = time_stepper_pt->time();
103 
104  Vector<double> deriv = derivative_function(t, u);
105  for(unsigned j=0, nj=deriv.size(); j<nj; j++)
106  {
107  // Get dudt approximation from time stepper
108  double dudt = time_stepper_pt->time_derivative(1, dat_pt, j);
109 
110  // Residual is difference between the exact derivative and our
111  // time stepper's derivative estimate.
112  residuals[j] = deriv[j] - dudt;
113  }
114  }
115 
117  DenseMatrix<double>& jacobian)
118  {
119  // Get residuals
121 
123  {
124  // get df/du jacobian
125  double t = internal_data_pt(0)->time_stepper_pt()->time();
126  Vector<double> dummy, u(nvalue(), 0.0);
127  internal_data_pt(0)->value(u);
128  Exact_solution_pt->jacobian(t, dummy, u, jacobian);
129 
130  // We need jacobian of residual = f - dudt so subtract diagonal
131  // (dudt)/du term.
132  const double a = internal_data_pt(0)->time_stepper_pt()->weight(1,0);
133  const unsigned n = nvalue();
134  for(unsigned i=0; i<n; i++)
135  {
136  jacobian(i, i) -= a;
137  }
138  }
139  else
140  {
141  // Use FD for jacobian
143  (residuals, jacobian, true);
144  }
145  }
146 
149  {
151  for(unsigned j=0, nj=nvalue(); j<nj; j++)
152  {
153  mm(j, j) = 1;
154  }
155  }
156 
157  /// Exact solution
158  Vector<double> exact_solution(const double& t) const
159  {
160 #ifdef PARANOID
161  if(Exact_solution_pt == 0)
162  {
163  throw OomphLibError("No exact solution function",
164  OOMPH_EXCEPTION_LOCATION,
165  OOMPH_CURRENT_FUNCTION);
166  }
167 #endif
168  Vector<double> dummy_x;
169  return (*Exact_solution_pt)(t, dummy_x);
170  }
171 
172  /// Exact solution
174  const Vector<double>& u)
175  {
176 #ifdef PARANOID
177  if(Exact_solution_pt == 0)
178  {
179  throw OomphLibError("No derivative function",
180  OOMPH_EXCEPTION_LOCATION,
181  OOMPH_CURRENT_FUNCTION);
182  }
183 #endif
184  Vector<double> dummy_x;
185  return Exact_solution_pt->derivative(t, dummy_x, u);
186  }
187 
189 
191  };
192 
193 }
194 
195 #endif
virtual void fill_in_contribution_to_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &mm)
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: ode_elements.h:147
A Generalised Element class.
Definition: elements.h:76
ODEElement()
Default constructor: null any pointers.
Definition: ode_elements.h:51
void fill_in_jacobian_from_internal_by_fd(Vector< double > &residuals, DenseMatrix< double > &jacobian, const bool &fd_all_data=false)
Calculate the contributions to the jacobian from the internal degrees of freedom using finite differe...
Definition: elements.cc:1064
Vector< double > exact_solution(const double &t) const
Exact solution.
Definition: ode_elements.h:158
cstr elem_len * i
Definition: cfortran.h:607
Element for integrating an initial value ODE.
Definition: ode_elements.h:45
double & time()
Return current value of continous time.
Definition: timesteppers.h:324
void build(TimeStepper *time_stepper_pt, SolutionFunctorBase *exact_solution_pt)
Store pointers, create internal data.
Definition: ode_elements.h:67
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
virtual Vector< double > derivative(const double &t, const Vector< double > &x, const Vector< double > &u) const =0
Call the derivative function.
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:448
unsigned nvalue() const
Definition: ode_elements.h:83
virtual ~ODEElement()
Definition: ode_elements.h:81
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:623
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:89
Vector< double > derivative_function(const double &t, const Vector< double > &u)
Exact solution.
Definition: ode_elements.h:173
virtual void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Add the elemental contribution to the jacobian matrix. and the residuals vector. Note that this funct...
Definition: ode_elements.h:116
virtual bool have_jacobian() const
Is a jacobian function implemented?
ODEElement(TimeStepper *time_stepper_pt, SolutionFunctorBase *exact_solution_pt)
Definition: ode_elements.h:60
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:246
unsigned add_internal_data(Data *const &data_pt, const bool &fd=true)
Add a (pointer to an) internal data object to the element and return the index required to obtain it ...
Definition: elements.cc:66
SolutionFunctorBase * Exact_solution_pt
Definition: ode_elements.h:188
virtual void jacobian(const double &t, const Vector< double > &x, const Vector< double > &u, DenseMatrix< double > &jacobian) const
void time_derivative(const unsigned &i, Data *const &data_pt, Vector< double > &deriv)
Evaluate i-th derivative of all values in Data and return in Vector deriv[].
Definition: timesteppers.h:474
double value(const unsigned &i) const
Return i-th stored value. This function is not virtual so that it can be inlined. This means that if ...
Definition: nodes.h:291
virtual void fill_in_contribution_to_residuals(Vector< double > &residuals)
Get residuals.
Definition: ode_elements.h:89
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition: timesteppers.h:219