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