gen_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 GeneralisedAdvection Diffusion elements
28
29namespace oomph
30{
31 /// 2D GeneralisedAdvection Diffusion elements
32
33
34 /// Default value for Peclet number
35 template<unsigned DIM>
37 0.0;
38
39 //======================================================================
40 /// Compute element residual Vector and/or element Jacobian matrix
41 ///
42 /// flag=1: compute both
43 /// flag=0: compute only residual Vector
44 ///
45 /// Pure version without hanging nodes
46 //======================================================================
47 template<unsigned DIM>
50 Vector<double>& residuals,
51 DenseMatrix<double>& jacobian,
52 DenseMatrix<double>& mass_matrix,
53 unsigned flag)
54 {
55 // Find out how many nodes there are
56 const unsigned n_node = nnode();
57
58 // Get the nodal index at which the unknown is stored
59 const unsigned u_nodal_index = u_index_cons_adv_diff();
60
61 // Set up memory for the shape and test functions
62 Shape psi(n_node), test(n_node);
63 DShape dpsidx(n_node, DIM), dtestdx(n_node, DIM);
64
65 // Set the value of n_intpt
66 const unsigned n_intpt = integral_pt()->nweight();
67
68 // Set the Vector to hold local coordinates
69 Vector<double> s(DIM);
70
71 // Get Peclet number
72 const double peclet = pe();
73
74 // Get the Peclet*Strouhal number
75 const double peclet_st = pe_st();
76
77 // Integers used to store the local equation number and local unknown
78 // indices for the residuals and jacobians
79 int local_eqn = 0, local_unknown = 0;
80
81 // Loop over the integration points
82 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
83 {
84 // Assign values of s
85 for (unsigned i = 0; i < DIM; i++) s[i] = integral_pt()->knot(ipt, i);
86
87 // Get the integral weight
88 double w = integral_pt()->weight(ipt);
89
90 // Call the derivatives of the shape and test functions
91 double J = dshape_and_dtest_eulerian_at_knot_cons_adv_diff(
92 ipt, psi, dpsidx, test, dtestdx);
93
94 // Premultiply the weights and the Jacobian
95 double W = w * J;
96
97 // Calculate local values of the solution and its derivatives
98 // Allocate
99 double interpolated_u = 0.0;
100 double dudt = 0.0;
101 Vector<double> interpolated_x(DIM, 0.0);
102 Vector<double> interpolated_dudx(DIM, 0.0);
103 Vector<double> mesh_velocity(DIM, 0.0);
104
105
106 // Calculate function value and derivatives:
107 //-----------------------------------------
108 // Loop over nodes
109 for (unsigned l = 0; l < n_node; l++)
110 {
111 // Get the value at the node
112 double u_value = raw_nodal_value(l, u_nodal_index);
113 interpolated_u += u_value * psi(l);
114 dudt += du_dt_cons_adv_diff(l) * psi(l);
115 // Loop over directions
116 for (unsigned j = 0; j < DIM; j++)
117 {
118 interpolated_x[j] += raw_nodal_position(l, j) * psi(l);
119 interpolated_dudx[j] += u_value * dpsidx(l, j);
120 }
121 }
122
123 // Mesh velocity?
124 if (!ALE_is_disabled)
125 {
126 for (unsigned l = 0; l < n_node; l++)
127 {
128 for (unsigned j = 0; j < DIM; j++)
129 {
130 mesh_velocity[j] += raw_dnodal_position_dt(l, j) * psi(l);
131 }
132 }
133 }
134
135
136 // Get source function
137 //-------------------
138 double source;
139 get_source_cons_adv_diff(ipt, interpolated_x, source);
140
141
142 // Get wind
143 //--------
144 Vector<double> wind(DIM);
145 get_wind_cons_adv_diff(ipt, s, interpolated_x, wind);
146
147 // Get the conserved wind (non-divergence free)
148 Vector<double> conserved_wind(DIM);
149 get_conserved_wind_cons_adv_diff(ipt, s, interpolated_x, conserved_wind);
150
151
152 // Get diffusivity tensor
153 DenseMatrix<double> D(DIM, DIM);
154 get_diff_cons_adv_diff(ipt, s, interpolated_x, D);
155
156 // Assemble residuals and Jacobian
157 //--------------------------------
158
159 // Loop over the test functions
160 for (unsigned l = 0; l < n_node; l++)
161 {
162 // Set the local equation number
163 local_eqn = nodal_local_eqn(l, u_nodal_index);
164
165 /*IF it's not a boundary condition*/
166 if (local_eqn >= 0)
167 {
168 // Add body force/source term and time derivative
169 residuals[local_eqn] -= (peclet_st * dudt + source) * test(l) * W;
170
171 // The Generalised Advection Diffusion bit itself
172 for (unsigned k = 0; k < DIM; k++)
173 {
174 // Terms that multiply the test function
175 // divergence-free wind
176 double tmp = peclet * wind[k];
177 // If the mesh is moving need to subtract the mesh velocity
178 if (!ALE_is_disabled)
179 {
180 tmp -= peclet_st * mesh_velocity[k];
181 }
182 tmp *= interpolated_dudx[k];
183
184 // Terms that multiply the derivative of the test function
185 // Advective term
186 double tmp2 = -conserved_wind[k] * interpolated_u;
187 // Now the diuffusive term
188 for (unsigned j = 0; j < DIM; j++)
189 {
190 tmp2 += D(k, j) * interpolated_dudx[j];
191 }
192 // Now construct the contribution to the residuals
193 residuals[local_eqn] -= (tmp * test(l) + tmp2 * dtestdx(l, k)) * W;
194 }
195
196 // Calculate the jacobian
197 //-----------------------
198 if (flag)
199 {
200 // Loop over the velocity shape functions again
201 for (unsigned l2 = 0; l2 < n_node; l2++)
202 {
203 // Set the number of the unknown
204 local_unknown = nodal_local_eqn(l2, u_nodal_index);
205
206 // If at a non-zero degree of freedom add in the entry
207 if (local_unknown >= 0)
208 {
209 // Mass matrix term
210 jacobian(local_eqn, local_unknown) -=
211 peclet_st * test(l) * psi(l2) *
212 node_pt(l2)->time_stepper_pt()->weight(1, 0) * W;
213
214 // Add the mass matrix term
215 if (flag == 2)
216 {
217 mass_matrix(local_eqn, local_unknown) +=
218 peclet_st * test(l) * psi(l2) * W;
219 }
220
221 // Add contribution to Elemental Matrix
222 for (unsigned k = 0; k < DIM; k++)
223 {
224 // Temporary term used in assembly
225 double tmp = -peclet * wind[k];
226 if (!ALE_is_disabled)
227 {
228 tmp -= peclet_st * mesh_velocity[k];
229 }
230 tmp *= dpsidx(l2, k);
231
232 double tmp2 = -conserved_wind[k] * psi(l2);
233 // Now the diffusive term
234 for (unsigned j = 0; j < DIM; j++)
235 {
236 tmp2 += D(k, j) * dpsidx(l2, j);
237 }
238
239 // Now assemble Jacobian term
240 jacobian(local_eqn, local_unknown) -=
241 (tmp * test(l) + tmp2 * dtestdx(l, k)) * W;
242 }
243 }
244 }
245 }
246 }
247 }
248
249
250 } // End of loop over integration points
251 }
252
253
254 //======================================================================
255 /// Self-test: Return 0 for OK
256 //======================================================================
257 template<unsigned DIM>
259 {
260 bool passed = true;
261
262 // Check lower-level stuff
263 if (FiniteElement::self_test() != 0)
264 {
265 passed = false;
266 }
267
268 // Return verdict
269 if (passed)
270 {
271 return 0;
272 }
273 else
274 {
275 return 1;
276 }
277 }
278
279
280 //======================================================================
281 /// Output function:
282 ///
283 /// x,y,u,w_x,w_y or x,y,z,u,w_x,w_y,w_z
284 ///
285 /// nplot points in each coordinate direction
286 //======================================================================
287 template<unsigned DIM>
289 std::ostream& outfile, const unsigned& nplot)
290 {
291 // Vector of local coordinates
292 Vector<double> s(DIM);
293
294
295 // Tecplot header info
296 outfile << tecplot_zone_string(nplot);
297
298 const unsigned n_node = this->nnode();
299 Shape psi(n_node);
300 DShape dpsidx(n_node, DIM);
301
302 // Loop over plot points
303 unsigned num_plot_points = nplot_points(nplot);
304 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
305 {
306 // Get local coordinates of plot point
307 get_s_plot(iplot, nplot, s);
308
309 // Get Eulerian coordinate of plot point
310 Vector<double> x(DIM);
311 interpolated_x(s, x);
312
313 for (unsigned i = 0; i < DIM; i++)
314 {
315 outfile << x[i] << " ";
316 }
317 outfile << interpolated_u_cons_adv_diff(s) << " ";
318
319 // Get the gradients
320 (void)this->dshape_eulerian(s, psi, dpsidx);
321 Vector<double> interpolated_dudx(DIM, 0.0);
322 for (unsigned n = 0; n < n_node; n++)
323 {
324 const double u_ = this->nodal_value(n, 0);
325 for (unsigned i = 0; i < DIM; i++)
326 {
327 interpolated_dudx[i] += u_ * dpsidx(n, i);
328 }
329 }
330
331 for (unsigned i = 0; i < DIM; i++)
332 {
333 outfile << interpolated_dudx[i] << " ";
334 }
335
336 // Get the wind
337 Vector<double> wind(DIM);
338 // Dummy integration point variable needed
339 unsigned ipt = 0;
340 get_wind_cons_adv_diff(ipt, s, x, wind);
341 for (unsigned i = 0; i < DIM; i++)
342 {
343 outfile << wind[i] << " ";
344 }
345 outfile << std::endl;
346 }
347
348 // Write tecplot footer (e.g. FE connectivity lists)
349 write_tecplot_zone_footer(outfile, nplot);
350 }
351
352
353 //======================================================================
354 /// C-style output function:
355 ///
356 /// x,y,u or x,y,z,u
357 ///
358 /// nplot points in each coordinate direction
359 //======================================================================
360 template<unsigned DIM>
362 FILE* file_pt, const unsigned& nplot)
363 {
364 // Vector of local coordinates
365 Vector<double> s(DIM);
366
367 // Tecplot header info
368 fprintf(file_pt, "%s", tecplot_zone_string(nplot).c_str());
369
370 // Loop over plot points
371 unsigned num_plot_points = nplot_points(nplot);
372 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
373 {
374 // Get local coordinates of plot point
375 get_s_plot(iplot, nplot, s);
376
377 for (unsigned i = 0; i < DIM; i++)
378 {
379 fprintf(file_pt, "%g ", interpolated_x(s, i));
380 }
381 fprintf(file_pt, "%g \n", interpolated_u_cons_adv_diff(s));
382 }
383
384 // Write tecplot footer (e.g. FE connectivity lists)
385 write_tecplot_zone_footer(file_pt, nplot);
386 }
387
388
389 //======================================================================
390 /// Output exact solution
391 ///
392 /// Solution is provided via function pointer.
393 /// Plot at a given number of plot points.
394 ///
395 /// x,y,u_exact or x,y,z,u_exact
396 //======================================================================
397 template<unsigned DIM>
399 std::ostream& outfile,
400 const unsigned& nplot,
402 {
403 // Vector of local coordinates
404 Vector<double> s(DIM);
405
406 // Vector for coordintes
407 Vector<double> x(DIM);
408
409 // Tecplot header info
410 outfile << tecplot_zone_string(nplot);
411
412 // Exact solution Vector (here a scalar)
413 Vector<double> exact_soln(1);
414
415 // Loop over plot points
416 unsigned num_plot_points = nplot_points(nplot);
417 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
418 {
419 // Get local coordinates of plot point
420 get_s_plot(iplot, nplot, s);
421
422 // Get x position as Vector
423 interpolated_x(s, x);
424
425 // Get exact solution at this point
426 (*exact_soln_pt)(x, exact_soln);
427
428 // Output x,y,...,u_exact
429 for (unsigned i = 0; i < DIM; i++)
430 {
431 outfile << x[i] << " ";
432 }
433 outfile << exact_soln[0] << std::endl;
434 }
435
436 // Write tecplot footer (e.g. FE connectivity lists)
437 write_tecplot_zone_footer(outfile, nplot);
438 }
439
440
441 //======================================================================
442 /// Validate against exact solution
443 ///
444 /// Solution is provided via function pointer.
445 /// Plot error at a given number of plot points.
446 ///
447 //======================================================================
448 template<unsigned DIM>
450 std::ostream& outfile,
452 double& error,
453 double& norm)
454 {
455 // Initialise
456 error = 0.0;
457 norm = 0.0;
458
459 // Vector of local coordinates
460 Vector<double> s(DIM);
461
462 // Vector for coordintes
463 Vector<double> x(DIM);
464
465 // Find out how many nodes there are in the element
466 unsigned n_node = nnode();
467
468 Shape psi(n_node);
469
470 // Set the value of n_intpt
471 unsigned n_intpt = integral_pt()->nweight();
472
473 // Tecplot header info
474 outfile << "ZONE" << std::endl;
475
476 // Exact solution Vector (here a scalar)
477 Vector<double> exact_soln(1);
478
479 // Loop over the integration points
480 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
481 {
482 // Assign values of s
483 for (unsigned i = 0; i < DIM; i++)
484 {
485 s[i] = integral_pt()->knot(ipt, i);
486 }
487
488 // Get the integral weight
489 double w = integral_pt()->weight(ipt);
490
491 // Get jacobian of mapping
492 double J = J_eulerian(s);
493
494 // Premultiply the weights and the Jacobian
495 double W = w * J;
496
497 // Get x position as Vector
498 interpolated_x(s, x);
499
500 // Get FE function value
501 double u_fe = interpolated_u_cons_adv_diff(s);
502
503 // Get exact solution at this point
504 (*exact_soln_pt)(x, exact_soln);
505
506 // Output x,y,...,error
507 for (unsigned i = 0; i < DIM; i++)
508 {
509 outfile << x[i] << " ";
510 }
511 outfile << exact_soln[0] << " " << exact_soln[0] - u_fe << std::endl;
512
513 // Add to error and norm
514 norm += exact_soln[0] * exact_soln[0] * W;
515 error += (exact_soln[0] - u_fe) * (exact_soln[0] - u_fe) * W;
516 }
517 }
518
519 //======================================================================
520 /// Calculate the integrated value of the unknown over the element
521 ///
522 //======================================================================
523 template<unsigned DIM>
525 {
526 // Initialise
527 double sum = 0.0;
528
529 // Vector of local coordinates
530 Vector<double> s(DIM);
531
532 // Find out how many nodes there are in the element
533 const unsigned n_node = nnode();
534
535 // Find the index at which the concentration is stored
536 const unsigned u_nodal_index = this->u_index_cons_adv_diff();
537
538 // Allocate memory for the shape functions
539 Shape psi(n_node);
540
541 // Set the value of n_intpt
542 const unsigned n_intpt = integral_pt()->nweight();
543
544 // Loop over the integration points
545 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
546 {
547 // Get the integral weight
548 const double w = integral_pt()->weight(ipt);
549
550 // Get the shape functions
551 this->shape_at_knot(ipt, psi);
552
553 // Calculate the concentration
554 double interpolated_u = 0.0;
555 for (unsigned l = 0; l < n_node; l++)
556 {
557 interpolated_u += this->nodal_value(l, u_nodal_index) * psi(l);
558 }
559
560 // Get jacobian of mapping
561 const double J = J_eulerian_at_knot(ipt);
562
563 // Add the values to the sum
564 sum += interpolated_u * w * J;
565 }
566
567 // return the sum
568 return sum;
569 }
570
571
572 //======================================================================
573 // Set the data for the number of Variables at each node
574 //======================================================================
575 template<unsigned DIM, unsigned NNODE_1D>
576 const unsigned
578
579 //====================================================================
580 // Force build of templates
581 //====================================================================
585
589
593
597
598
599} // 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(* 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 all elements that solve the Advection Diffusion equations in conservative form using isop...
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
static double Default_peclet_number
Static default value for the Peclet number.
void output(std::ostream &outfile)
Output with default number of plot points.
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.
double integrate_u()
Integrate the concentration over the element.
virtual void fill_in_generic_residual_contribution_cons_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.
//////////////////////////////////////////////////////////////////////// ////////////////////////////...
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
//////////////////////////////////////////////////////////////////// ////////////////////////////////...