advection_diffusion_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 functions for Advection Diffusion elements
28
29namespace oomph
30{
31 /// 2D Advection Diffusion elements
32
33
34 /// Default value for Peclet number
35 template<unsigned DIM>
37
38 //======================================================================
39 /// Compute element residual Vector and/or element Jacobian matrix
40 ///
41 /// flag=1: compute both
42 /// flag=0: compute only residual Vector
43 ///
44 /// Pure version without hanging nodes
45 //======================================================================
46 template<unsigned DIM>
49 Vector<double>& residuals,
50 DenseMatrix<double>& jacobian,
51 DenseMatrix<double>& mass_matrix,
52 unsigned flag)
53 {
54 // Find out how many nodes there are
55 unsigned n_node = nnode();
56
57 // Get the nodal index at which the unknown is stored
58 unsigned u_nodal_index = u_index_adv_diff();
59
60 // Set up memory for the shape and test functions
61 Shape psi(n_node), test(n_node);
62 DShape dpsidx(n_node, DIM), dtestdx(n_node, DIM);
63
64 // Set the value of n_intpt
65 unsigned n_intpt = integral_pt()->nweight();
66
67 // Set the Vector to hold local coordinates
68 Vector<double> s(DIM);
69
70 // Get Peclet number
71 double peclet = pe();
72
73 // Get the Peclet*Strouhal number
74 double peclet_st = pe_st();
75
76 // Integers used to store the local equation number and local unknown
77 // indices for the residuals and jacobians
78 int local_eqn = 0, local_unknown = 0;
79
80 // Loop over the integration points
81 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
82 {
83 // Assign values of s
84 for (unsigned i = 0; i < DIM; i++) s[i] = integral_pt()->knot(ipt, i);
85
86 // Get the integral weight
87 double w = integral_pt()->weight(ipt);
88
89 // Call the derivatives of the shape and test functions
90 double J = dshape_and_dtest_eulerian_at_knot_adv_diff(
91 ipt, psi, dpsidx, test, dtestdx);
92
93 // Premultiply the weights and the Jacobian
94 double W = w * J;
95
96 // Calculate local values of the solution and its derivatives
97 // Allocate
98 double interpolated_u = 0.0;
99 double dudt = 0.0;
100 Vector<double> interpolated_x(DIM, 0.0);
101 Vector<double> interpolated_dudx(DIM, 0.0);
102 Vector<double> mesh_velocity(DIM, 0.0);
103
104
105 // Calculate function value and derivatives:
106 //-----------------------------------------
107 // Loop over nodes
108 for (unsigned l = 0; l < n_node; l++)
109 {
110 // Get the value at the node
111 double u_value = raw_nodal_value(l, u_nodal_index);
112 interpolated_u += u_value * psi(l);
113 dudt += du_dt_adv_diff(l) * psi(l);
114 // Loop over directions
115 for (unsigned j = 0; j < DIM; j++)
116 {
117 interpolated_x[j] += raw_nodal_position(l, j) * psi(l);
118 interpolated_dudx[j] += u_value * dpsidx(l, j);
119 }
120 }
121
122 // Mesh velocity?
123 if (!ALE_is_disabled)
124 {
125 for (unsigned l = 0; l < n_node; l++)
126 {
127 for (unsigned j = 0; j < DIM; j++)
128 {
129 mesh_velocity[j] += raw_dnodal_position_dt(l, j) * psi(l);
130 }
131 }
132 }
133
134
135 // Get source function
136 //-------------------
137 double source;
138 get_source_adv_diff(ipt, interpolated_x, source);
139
140
141 // Get wind
142 //--------
143 Vector<double> wind(DIM);
144 get_wind_adv_diff(ipt, s, interpolated_x, wind);
145
146 // Assemble residuals and Jacobian
147 //--------------------------------
148
149 // Loop over the test functions
150 for (unsigned l = 0; l < n_node; l++)
151 {
152 // Set the local equation number
153 local_eqn = nodal_local_eqn(l, u_nodal_index);
154
155 /*IF it's not a boundary condition*/
156 if (local_eqn >= 0)
157 {
158 // Add body force/source term and time derivative
159 residuals[local_eqn] -= (peclet_st * dudt + source) * test(l) * W;
160
161 // The Advection Diffusion bit itself
162 for (unsigned k = 0; k < DIM; k++)
163 {
164 // Terms that multiply the test function
165 double tmp = peclet * wind[k];
166 // If the mesh is moving need to subtract the mesh velocity
167 if (!ALE_is_disabled)
168 {
169 tmp -= peclet_st * mesh_velocity[k];
170 }
171 // Now construct the contribution to the residuals
172 residuals[local_eqn] -=
173 interpolated_dudx[k] * (tmp * test(l) + dtestdx(l, k)) * W;
174 }
175
176 // Calculate the jacobian
177 //-----------------------
178 if (flag)
179 {
180 // Loop over the velocity shape functions again
181 for (unsigned l2 = 0; l2 < n_node; l2++)
182 {
183 // Set the number of the unknown
184 local_unknown = nodal_local_eqn(l2, u_nodal_index);
185
186 // If at a non-zero degree of freedom add in the entry
187 if (local_unknown >= 0)
188 {
189 // Mass matrix term
190 jacobian(local_eqn, local_unknown) -=
191 peclet_st * test(l) * psi(l2) *
192 node_pt(l2)->time_stepper_pt()->weight(1, 0) * W;
193
194 // Add the mass matrix term
195 if (flag == 2)
196 {
197 mass_matrix(local_eqn, local_unknown) +=
198 peclet_st * test(l) * psi(l2) * W;
199 }
200
201 // Add contribution to Elemental Matrix
202 for (unsigned i = 0; i < DIM; i++)
203 {
204 // Temporary term used in assembly
205 double tmp = peclet * wind[i];
206 if (!ALE_is_disabled) tmp -= peclet_st * mesh_velocity[i];
207 // Now assemble Jacobian term
208 jacobian(local_eqn, local_unknown) -=
209 dpsidx(l2, i) * (tmp * test(l) + dtestdx(l, i)) * W;
210 }
211 }
212 }
213 }
214 }
215 }
216
217
218 } // End of loop over integration points
219 }
220
221
222 //======================================================================
223 /// Self-test: Return 0 for OK
224 //======================================================================
225 template<unsigned DIM>
227 {
228 bool passed = true;
229
230 // Check lower-level stuff
231 if (FiniteElement::self_test() != 0)
232 {
233 passed = false;
234 }
235
236 // Return verdict
237 if (passed)
238 {
239 return 0;
240 }
241 else
242 {
243 return 1;
244 }
245 }
246
247
248 //======================================================================
249 /// Output function:
250 ///
251 /// x,y,u,w_x,w_y or x,y,z,u,w_x,w_y,w_z
252 ///
253 /// nplot points in each coordinate direction
254 //======================================================================
255 template<unsigned DIM>
257 const unsigned& nplot)
258 {
259 // Vector of local coordinates
260 Vector<double> s(DIM);
261
262
263 // Tecplot header info
264 outfile << tecplot_zone_string(nplot);
265
266 // Loop over plot points
267 unsigned num_plot_points = nplot_points(nplot);
268 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
269 {
270 // Get local coordinates of plot point
271 get_s_plot(iplot, nplot, s);
272
273 // Get Eulerian coordinate of plot point
274 Vector<double> x(DIM);
275 interpolated_x(s, x);
276
277 for (unsigned i = 0; i < DIM; i++)
278 {
279 outfile << x[i] << " ";
280 }
281 outfile << interpolated_u_adv_diff(s) << " ";
282
283 // Get the wind
284 Vector<double> wind(DIM);
285 // Dummy ipt argument needed... ?
286 unsigned ipt = 0;
287 get_wind_adv_diff(ipt, s, x, wind);
288 for (unsigned i = 0; i < DIM; i++)
289 {
290 outfile << wind[i] << " ";
291 }
292 outfile << std::endl;
293 }
294
295 // Write tecplot footer (e.g. FE connectivity lists)
296 write_tecplot_zone_footer(outfile, nplot);
297 }
298
299
300 //======================================================================
301 /// C-style output function:
302 ///
303 /// x,y,u or x,y,z,u
304 ///
305 /// nplot points in each coordinate direction
306 //======================================================================
307 template<unsigned DIM>
309 const unsigned& nplot)
310 {
311 // Vector of local coordinates
312 Vector<double> s(DIM);
313
314 // Tecplot header info
315 fprintf(file_pt, "%s", tecplot_zone_string(nplot).c_str());
316
317 // Loop over plot points
318 unsigned num_plot_points = nplot_points(nplot);
319 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
320 {
321 // Get local coordinates of plot point
322 get_s_plot(iplot, nplot, s);
323
324 for (unsigned i = 0; i < DIM; i++)
325 {
326 fprintf(file_pt, "%g ", interpolated_x(s, i));
327 }
328 fprintf(file_pt, "%g \n", interpolated_u_adv_diff(s));
329 }
330
331 // Write tecplot footer (e.g. FE connectivity lists)
332 write_tecplot_zone_footer(file_pt, nplot);
333 }
334
335
336 //======================================================================
337 /// Output exact solution
338 ///
339 /// Solution is provided via function pointer.
340 /// Plot at a given number of plot points.
341 ///
342 /// x,y,u_exact or x,y,z,u_exact
343 //======================================================================
344 template<unsigned DIM>
346 std::ostream& outfile,
347 const unsigned& nplot,
349 {
350 // Vector of local coordinates
351 Vector<double> s(DIM);
352
353 // Vector for coordintes
354 Vector<double> x(DIM);
355
356 // Tecplot header info
357 outfile << tecplot_zone_string(nplot);
358
359 // Exact solution Vector (here a scalar)
360 Vector<double> exact_soln(1);
361
362 // Loop over plot points
363 unsigned num_plot_points = nplot_points(nplot);
364 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
365 {
366 // Get local coordinates of plot point
367 get_s_plot(iplot, nplot, s);
368
369 // Get x position as Vector
370 interpolated_x(s, x);
371
372 // Get exact solution at this point
373 (*exact_soln_pt)(x, exact_soln);
374
375 // Output x,y,...,u_exact
376 for (unsigned i = 0; i < DIM; i++)
377 {
378 outfile << x[i] << " ";
379 }
380 outfile << exact_soln[0] << std::endl;
381 }
382
383 // Write tecplot footer (e.g. FE connectivity lists)
384 write_tecplot_zone_footer(outfile, nplot);
385 }
386
387
388 //======================================================================
389 /// Validate against exact solution
390 ///
391 /// Solution is provided via function pointer.
392 /// Plot error at a given number of plot points.
393 ///
394 //======================================================================
395 template<unsigned DIM>
397 std::ostream& outfile,
399 double& error,
400 double& norm)
401 {
402 // Initialise
403 error = 0.0;
404 norm = 0.0;
405
406 // Vector of local coordinates
407 Vector<double> s(DIM);
408
409 // Vector for coordintes
410 Vector<double> x(DIM);
411
412 // Find out how many nodes there are in the element
413 unsigned n_node = nnode();
414
415 Shape psi(n_node);
416
417 // Set the value of n_intpt
418 unsigned n_intpt = integral_pt()->nweight();
419
420 // Tecplot header info
421 outfile << "ZONE" << std::endl;
422
423 // Exact solution Vector (here a scalar)
424 Vector<double> exact_soln(1);
425
426 // Loop over the integration points
427 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
428 {
429 // Assign values of s
430 for (unsigned i = 0; i < DIM; i++)
431 {
432 s[i] = integral_pt()->knot(ipt, i);
433 }
434
435 // Get the integral weight
436 double w = integral_pt()->weight(ipt);
437
438 // Get jacobian of mapping
439 double J = J_eulerian(s);
440
441 // Premultiply the weights and the Jacobian
442 double W = w * J;
443
444 // Get x position as Vector
445 interpolated_x(s, x);
446
447 // Get FE function value
448 double u_fe = interpolated_u_adv_diff(s);
449
450 // Get exact solution at this point
451 (*exact_soln_pt)(x, exact_soln);
452
453 // Output x,y,...,error
454 for (unsigned i = 0; i < DIM; i++)
455 {
456 outfile << x[i] << " ";
457 }
458 outfile << exact_soln[0] << " " << exact_soln[0] - u_fe << std::endl;
459
460 // Add to error and norm
461 norm += exact_soln[0] * exact_soln[0] * W;
462 error += (exact_soln[0] - u_fe) * (exact_soln[0] - u_fe) * W;
463 }
464 }
465
466
467 //======================================================================
468 // Set the data for the number of Variables at each node
469 //======================================================================
470 template<unsigned DIM, unsigned NNODE_1D>
472
473 //====================================================================
474 // Force build of templates
475 //====================================================================
479
483
487
488
489} // namespace oomph
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
A class for all elements that solve the Advection Diffusion equations using isoparametric elements.
void output(std::ostream &outfile)
Output with default number of plot points.
virtual void fill_in_generic_residual_contribution_adv_diff(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
Add the element's contribution to its residual vector only (if flag=and/or element Jacobian matrix.
void output_fct(std::ostream &outfile, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output exact soln: x,y,u_exact or x,y,z,u_exact at nplot^DIM plot points.
unsigned self_test()
Self-test: Return 0 for OK.
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as .
Definition: elements.h:1759
virtual unsigned self_test()
Self-test: Check inversion of element & do self-test for GeneralisedElement. Return 0 if OK.
Definition: elements.cc:4440
//////////////////////////////////////////////////////////////////////// ////////////////////////////...
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
//////////////////////////////////////////////////////////////////// ////////////////////////////////...