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// 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//====================================================================
27
28namespace oomph
29{
31 std::ostream& outfile,
32 const unsigned& n_plot)
33 {
34 // Find out how many nodes there are
35 const unsigned n_node = nnode();
36
37 // Set up memory for the shape and test functions
38 Shape psi(n_node);
39 DShape dpsidt(n_node, 1);
40
42
43 // Cast the element
44 PeriodicOrbitBaseElement* const base_el_pt =
45 dynamic_cast<PeriodicOrbitBaseElement*>(elem_pt);
46
47 const double inverse_timescale = this->omega();
48
49 // Loop over the plot point
50 for (unsigned i = 0; i < n_plot; i++)
51 {
52 // Get the coordinate
53 s[0] = -1.0 + 2.0 * i / ((double)(n_plot - 1));
54
55 (void)dshape_eulerian(s, psi, dpsidt);
56
57 // Sort out the timestepper weights and the time in here
58 // Calculate the time
59 double interpolated_time = 0.0;
60 for (unsigned l = 0; l < n_node; l++)
61 {
62 interpolated_time += raw_nodal_position(l, 0) * psi(l);
63 }
64
65
66 // Need to multiply by period and the set global time
67 this->time_pt()->time() = interpolated_time;
68
69 // Set the weights of the timestepper
70 this->set_timestepper_weights(psi, dpsidt);
71
72 base_el_pt->spacetime_output(
73 outfile, n_plot, interpolated_time / inverse_timescale);
74 }
75 }
76
78 PeriodicOrbitAssemblyHandlerBase* const& assembly_handler_pt,
79 GeneralisedElement* const& elem_pt,
80 Vector<double>& residuals,
81 DenseMatrix<double>& jacobian,
82 const unsigned& flag)
83 {
84 // A simple integration loop
85 // Find out how many nodes there are
86 const unsigned n_node = nnode();
87
88 // Set up memory for the shape and test functions
89 Shape psi(n_node), test(n_node);
90 DShape dpsidt(n_node, 1), dtestdt(n_node, 1);
91
92 // Set the value of n_intpt
93 const unsigned n_intpt = integral_pt()->nweight();
94
95 // Integers to store the local equation and unknown numbers
96 int local_eqn = 0, local_unknown = 0;
97
98 // Storage for the underlying element's residuals
99 const unsigned n_elem_dof = elem_pt->ndof();
100 Vector<double> el_residuals(n_elem_dof);
101 DenseMatrix<double> el_mass;
102 DenseMatrix<double> el_jacobian;
103
104 if (flag)
105 {
106 el_mass.resize(n_elem_dof, n_elem_dof);
107 el_jacobian.resize(n_elem_dof, n_elem_dof);
108 }
109
110 // Storage for the current value of the unkowns
111 Vector<double> u(n_elem_dof);
112 // and the derivatives
113 Vector<double> du_dt(n_elem_dof);
114 // And the previous derivative
115 Vector<double> du_dt_old(n_elem_dof);
116 // Storage for the inner product matrix
117 DenseMatrix<double> inner_product(n_elem_dof);
118 // Cast the element
119 PeriodicOrbitBaseElement* const base_el_pt =
120 dynamic_cast<PeriodicOrbitBaseElement*>(elem_pt);
121 // Get the inner product matrix
122 base_el_pt->get_inner_product_matrix(inner_product);
123
124 // Storage for "all" unknowns
125 Vector<double> all_current_unknowns;
126 Vector<double> all_previous_unknowns;
127
128 // Get the value of omega
129 const double inverse_timescale = this->omega();
130
131 // Loop over the integration points
132 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
133 {
134 // Get the integral weight
135 double w = integral_pt()->weight(ipt);
136
137 // Call the derivatives of the shape and test functions
139 ipt, psi, dpsidt, test, dtestdt);
140
141 // Premultiply the weights and the Jacobian
142 double W = w * J;
143
144 // Sort out the timestepper weights and the time in here
145 // Calculate the time
146 double interpolated_time = 0.0;
147 for (unsigned l = 0; l < n_node; l++)
148 {
149 interpolated_time += raw_nodal_position(l, 0) * psi(l);
150 }
151
152
153 // Need to multiply by period and the set global time
154 this->time_pt()->time() = interpolated_time / inverse_timescale;
155 // Set the weights of the timestepper
156 this->set_timestepper_weights(psi, dpsidt);
157
158 // Get the residuals of the element or residuals and jacobian
159 if (flag)
160 {
162 el_residuals, el_jacobian, el_mass);
163 }
164 else
165 {
166 elem_pt->get_residuals(el_residuals);
167 }
168
169 // Get the current value of the unknown time derivatives
170 base_el_pt->get_non_external_ddofs_dt(du_dt);
171
172 // Multiply the elemental residuals by the appropriate test functions
173 for (unsigned l = 0; l < n_node; l++)
174 {
175 // Get the local equation number
176 local_eqn = nodal_local_eqn(l, 0);
177 // If it's not a boundary condition (it will never be)
178 if (local_eqn >= 0)
179 {
180 // Work out the offset which is the GLOBAL time unknown
181 // multiplied by the number of elemental unknowns
182 unsigned offset = this->eqn_number(local_eqn) * n_elem_dof;
183 // Add to the appropriate residuals
184 for (unsigned i = 0; i < n_elem_dof; i++)
185 {
186 residuals[i + offset] += el_residuals[i] * psi(l) * W;
187 }
188
189
190 // Add the jacobian contributions
191 if (flag)
192 {
193 // The form of the equations is -M du/dt + R = 0
194
195 // Now add the contribution to the jacobian
196 for (unsigned l2 = 0; l2 < n_node; l2++)
197 {
198 local_unknown = nodal_local_eqn(l2, 0);
199 if (local_unknown >= 0)
200 {
201 // Work out the second offset
202 unsigned offset2 = this->eqn_number(local_unknown) * n_elem_dof;
203 // Add to the appropriate jacobian terms
204 for (unsigned i = 0; i < n_elem_dof; i++)
205 {
206 for (unsigned j = 0; j < n_elem_dof; j++)
207 {
208 // Add in the Jacobian terms
209 jacobian(i + offset, j + offset2) +=
210 el_jacobian(i, j) * psi(l2) * psi(l) * W;
211
212 // Add the time derivative terms,
213 jacobian(i + offset, j + offset2) -=
214 el_mass(i, j) * dpsidt(l2, 0) * inverse_timescale *
215 psi(l) * W;
216 }
217 }
218 }
219 }
220
221
222 // Add the variation of the period
223 for (unsigned i = 0; i < n_elem_dof; i++)
224 {
225 for (unsigned j = 0; j < n_elem_dof; j++)
226 {
227 jacobian(i + offset, Ntstorage * n_elem_dof) -=
228 el_mass(i, j) * (du_dt[j] / inverse_timescale) * psi(l) * W;
229 }
230 }
231 } // End of Jacobian flag
232 }
233 }
234
235 // Sort out the phase condition
236
237 // Get the current value of the unknowns
238 base_el_pt->get_non_external_dofs(u);
239
240 // Now get the unknowns required by the assembly handler
241 // i.e. including all values throughout the period for backup
242 assembly_handler_pt->get_dofs_for_element(elem_pt, all_current_unknowns);
243 // Get the previous values as stored by the assembly handler
244 assembly_handler_pt->get_previous_dofs_for_element(elem_pt,
245 all_previous_unknowns);
246
247 // Now set the elemental values to the previous values
248 assembly_handler_pt->set_dofs_for_element(elem_pt, all_previous_unknowns);
249 // Get the previous time derivatives
250 base_el_pt->get_non_external_ddofs_dt(du_dt_old);
251 // Reset the element's values to the current
252 assembly_handler_pt->set_dofs_for_element(elem_pt, all_current_unknowns);
253
254
255 // Assemble the inner product
256 double sum = 0.0;
257 for (unsigned i = 0; i < n_elem_dof; i++)
258 {
259 for (unsigned j = 0; j < n_elem_dof; j++)
260 {
261 sum += u[i] * inner_product(i, j) * du_dt_old[j];
262 }
263 }
264
265 // Add to the residuals
266 residuals[Ntstorage * n_elem_dof] += sum * W;
267
268 // Sort out the jacobian
269 if (flag)
270 {
271 // Loop over the unknown time points
272 for (unsigned l2 = 0; l2 < n_node; l2++)
273 {
274 // Get the local unknown
275 local_unknown = nodal_local_eqn(l2, 0);
276 if (local_unknown >= 0)
277 {
278 // Work out the offset
279 unsigned offset2 = this->eqn_number(local_unknown) * n_elem_dof;
280 // Now add in the appropriate jacobian terms
281 for (unsigned i2 = 0; i2 < n_elem_dof; i2++)
282 {
283 double sum2 = 0.0;
284 for (unsigned j = 0; j < n_elem_dof; j++)
285 {
286 sum2 += inner_product(i2, j) * du_dt_old[j];
287 }
288 jacobian(Ntstorage * n_elem_dof, i2 + offset2) +=
289 psi(l2) * sum2 * W;
290 }
291 }
292 }
293 } // End of jacobian calculation
294
295 } // End of loop over time period integration points
296 }
297
298
299} // namespace oomph
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
void resize(const unsigned long &n)
Resize to a square nxn matrix; any values already present will be transfered.
Definition: matrices.h:498
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
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 nnode() const
Return the number of nodes.
Definition: elements.h:2210
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:3298
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....
Definition: elements.cc:1686
A Generalised Element class.
Definition: elements.h:73
unsigned ndof() const
Return the number of equations/dofs in the element.
Definition: elements.h:835
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:704
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:980
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:1016
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.
=============================================================== Base class to avoid template complica...
virtual void set_dofs_for_element(GeneralisedElement *const elem_pt, Vector< double > const &dofs)=0
virtual void get_dofs_for_element(GeneralisedElement *const elem_pt, Vector< double > &dofs)=0
virtual void get_previous_dofs_for_element(GeneralisedElement *const elem_pt, Vector< double > &dofs)=0
virtual void get_non_external_ddofs_dt(Vector< double > &du_dt)
Interface to get the current value of the time derivative of all (internal and shared) unknowns.
virtual void spacetime_output(std::ostream &outilfe, const unsigned &Nplot, const double &time=0.0)
virtual void get_non_external_dofs(Vector< double > &u)
Interface to get the current value of all (internal and shared) unknowns.
virtual void get_inner_product_matrix(DenseMatrix< double > &inner_product)
Get the inner product matrix.
void orbit_output(GeneralisedElement *const &elem_pt, std::ostream &outfile, const unsigned &n_plot)
Time *& time_pt()
Retun the pointer to the global time.
void set_timestepper_weights(const Shape &psi, const DShape &dpsidt)
Set the timestepper weights.
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 ...
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!
double omega()
Return the frequency.
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
double & time()
Return the current value of the continuous time.
Definition: timesteppers.h:123
//////////////////////////////////////////////////////////////////// ////////////////////////////////...