gen_axisym_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 GeneralisedAxisymAdvection Diffusion elements
28
29namespace oomph
30{
31 /// 2D GeneralisedAdvection Diffusion elements
32
33
34 /// Default value for Peclet number
36 0.0;
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 //======================================================================
48 Vector<double>& residuals,
49 DenseMatrix<double>& jacobian,
50 DenseMatrix<double>& mass_matrix,
51 unsigned flag)
52 {
53 // Find out how many nodes there are
54 const unsigned n_node = nnode();
55
56 // Get the nodal index at which the unknown is stored
57 const unsigned u_nodal_index = u_index_cons_axisym_adv_diff();
58
59 // Set up memory for the shape and test functions
60 Shape psi(n_node), test(n_node);
61 DShape dpsidx(n_node, 2), dtestdx(n_node, 2);
62
63 // Set the value of n_intpt
64 const unsigned n_intpt = integral_pt()->nweight();
65
66 // Set the Vector to hold local coordinates
68
69 // Get Peclet number
70 const double peclet = pe();
71
72 // Get the Peclet*Strouhal number
73 const double peclet_st = pe_st();
74
75 // Integers used to store the local equation number and local unknown
76 // indices for the residuals and jacobians
77 int local_eqn = 0, local_unknown = 0;
78
79 // Loop over the integration points
80 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
81 {
82 // Assign values of s
83 for (unsigned i = 0; i < 2; i++) s[i] = integral_pt()->knot(ipt, i);
84
85 // Get the integral weight
86 double w = integral_pt()->weight(ipt);
87
88 // Call the derivatives of the shape and test functions
90 ipt, psi, dpsidx, test, dtestdx);
91
92 // Premultiply the weights and the Jacobian
93 double W = w * J;
94
95 // Calculate local values of the solution and its derivatives
96 // Allocate
97 double interpolated_u = 0.0;
98 double dudt = 0.0;
100 Vector<double> interpolated_dudx(2, 0.0);
101 Vector<double> mesh_velocity(2, 0.0);
102
103
104 // Calculate function value and derivatives:
105 //-----------------------------------------
106 // Loop over nodes
107 for (unsigned l = 0; l < n_node; l++)
108 {
109 // Get the value at the node
110 double u_value = raw_nodal_value(l, u_nodal_index);
111 interpolated_u += u_value * psi(l);
112 dudt += du_dt_cons_axisym_adv_diff(l) * psi(l);
113 // Loop over directions
114 for (unsigned j = 0; j < 2; j++)
115 {
116 interpolated_x[j] += raw_nodal_position(l, j) * psi(l);
117 interpolated_dudx[j] += u_value * dpsidx(l, j);
118 }
119 }
120
121 // Mesh velocity?
122 if (!ALE_is_disabled)
123 {
124 for (unsigned l = 0; l < n_node; l++)
125 {
126 for (unsigned j = 0; j < 2; j++)
127 {
128 mesh_velocity[j] += raw_dnodal_position_dt(l, j) * psi(l);
129 }
130 }
131 }
132
133
134 // Get source function
135 //-------------------
136 double source;
138
139
140 // Get wind (three components because this could come from a NS
141 // computation)
142 //--------
145
146 // Get the conserved wind (non-divergence free)
147 Vector<double> conserved_wind(3);
149 ipt, s, interpolated_x, conserved_wind);
150
151 // Get diffusivity tensor
154
155 // r is the first position component
156 double r = interpolated_x[0];
157
158 // Assemble residuals and Jacobian
159 //--------------------------------
160
161 // Loop over the test functions
162 for (unsigned l = 0; l < n_node; l++)
163 {
164 // Set the local equation number
165 local_eqn = nodal_local_eqn(l, u_nodal_index);
166
167 /*IF it's not a boundary condition*/
168 if (local_eqn >= 0)
169 {
170 // Add body force/source term and time derivative
171 residuals[local_eqn] -= (peclet_st * dudt + source) * r * test(l) * W;
172
173 // The Generalised Advection Diffusion bit itself
174 // Only loop over the non azimuthal directions
175 for (unsigned k = 0; k < 2; k++)
176 {
177 // Terms that multiply the test function
178 // divergence-free wind
179 double tmp = peclet * wind[k];
180 // If the mesh is moving need to subtract the mesh velocity
181 if (!ALE_is_disabled)
182 {
183 tmp -= peclet_st * mesh_velocity[k];
184 }
185 tmp *= interpolated_dudx[k];
186
187 // Terms that multiply the derivative of the test function
188 // Advective term
189 double tmp2 = -conserved_wind[k] * interpolated_u;
190 // Now the diuffusive term
191 for (unsigned j = 0; j < 2; j++)
192 {
193 tmp2 += D(k, j) * interpolated_dudx[j];
194 }
195 // Now construct the contribution to the residuals
196 residuals[local_eqn] -=
197 (tmp * test(l) + tmp2 * dtestdx(l, k)) * r * W;
198 }
199
200 // Calculate the jacobian
201 //-----------------------
202 if (flag)
203 {
204 // Loop over the velocity shape functions again
205 for (unsigned l2 = 0; l2 < n_node; l2++)
206 {
207 // Set the number of the unknown
208 local_unknown = nodal_local_eqn(l2, u_nodal_index);
209
210 // If at a non-zero degree of freedom add in the entry
211 if (local_unknown >= 0)
212 {
213 // Mass matrix term
214 jacobian(local_eqn, local_unknown) -=
215 peclet_st * test(l) * psi(l2) *
216 node_pt(l2)->time_stepper_pt()->weight(1, 0) * r * W;
217
218 // Add the mass matrix term
219 if (flag == 2)
220 {
221 mass_matrix(local_eqn, local_unknown) +=
222 peclet_st * test(l) * psi(l2) * r * W;
223 }
224
225 // Add contribution to Elemental Matrix
226 for (unsigned k = 0; k < 2; k++)
227 {
228 // Temporary term used in assembly
229 double tmp = -peclet * wind[k];
230 if (!ALE_is_disabled)
231 {
232 tmp -= peclet_st * mesh_velocity[k];
233 }
234 tmp *= dpsidx(l2, k);
235
236 double tmp2 = -conserved_wind[k] * psi(l2);
237 // Now the diffusive term
238 for (unsigned j = 0; j < 2; j++)
239 {
240 tmp2 += D(k, j) * dpsidx(l2, j);
241 }
242
243 // Now assemble Jacobian term
244 jacobian(local_eqn, local_unknown) -=
245 (tmp * test(l) + tmp2 * dtestdx(l, k)) * r * W;
246 }
247 }
248 }
249 }
250 }
251 }
252
253
254 } // End of loop over integration points
255 }
256
257
258 //======================================================================
259 /// Self-test: Return 0 for OK
260 //======================================================================
262 {
263 bool passed = true;
264
265 // Check lower-level stuff
266 if (FiniteElement::self_test() != 0)
267 {
268 passed = false;
269 }
270
271 // Return verdict
272 if (passed)
273 {
274 return 0;
275 }
276 else
277 {
278 return 1;
279 }
280 }
281
282
283 //======================================================================
284 /// Output function:
285 ///
286 /// r,z,u,w_r,w_z
287 ///
288 /// nplot points in each coordinate direction
289 //======================================================================
291 std::ostream& outfile, const unsigned& nplot)
292 {
293 // Vector of local coordinates
294 Vector<double> s(2);
295
296 // Tecplot header info
297 outfile << tecplot_zone_string(nplot);
298
299 const unsigned n_node = this->nnode();
300 Shape psi(n_node);
301 DShape dpsidx(n_node, 2);
302
303 // Loop over plot points
304 unsigned num_plot_points = nplot_points(nplot);
305 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
306 {
307 // Get local coordinates of plot point
308 get_s_plot(iplot, nplot, s);
309
310 // Get Eulerian coordinate of plot point
311 Vector<double> x(2);
312 interpolated_x(s, x);
313
314 for (unsigned i = 0; i < 2; i++)
315 {
316 outfile << x[i] << " ";
317 }
318 outfile << interpolated_u_cons_axisym_adv_diff(s) << " ";
319
320 // Get the gradients
321 /*(void)this->dshape_eulerian(s,psi,dpsidx);
322 Vector<double> interpolated_dudx(2,0.0);
323 for(unsigned n=0;n<n_node;n++)
324 {
325 const double u_ = this->nodal_value(n,0);
326 for(unsigned i=0;i<2;i++)
327 {
328 interpolated_dudx[i] += u_*dpsidx(n,i);
329 }
330 }
331
332 for(unsigned i=0;i<2;i++)
333 {
334 outfile << interpolated_dudx[i] << " ";
335 }*/
336
337 // Get the wind
339 // Dummy integration point variable needed
340 unsigned ipt = 0;
342 for (unsigned i = 0; i < 3; i++)
343 {
344 outfile << wind[i] << " ";
345 }
346 outfile << std::endl;
347 }
348
349 // Write tecplot footer (e.g. FE connectivity lists)
350 write_tecplot_zone_footer(outfile, nplot);
351 }
352
353
354 //======================================================================
355 /// C-style output function:
356 ///
357 /// r,z,u
358 ///
359 /// nplot points in each coordinate direction
360 //======================================================================
362 FILE* file_pt, const unsigned& nplot)
363 {
364 // Vector of local coordinates
365 Vector<double> s(2);
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 < 2; i++)
378 {
379 fprintf(file_pt, "%g ", interpolated_x(s, i));
380 }
381 fprintf(file_pt, "%g \n", interpolated_u_cons_axisym_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 /// r,z,,u_exact
396 //======================================================================
398 std::ostream& outfile,
399 const unsigned& nplot,
401 {
402 // Vector of local coordinates
403 Vector<double> s(2);
404
405 // Vector for coordintes
406 Vector<double> x(2);
407
408 // Tecplot header info
409 outfile << tecplot_zone_string(nplot);
410
411 // Exact solution Vector (here a scalar)
412 Vector<double> exact_soln(1);
413
414 // Loop over plot points
415 unsigned num_plot_points = nplot_points(nplot);
416 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
417 {
418 // Get local coordinates of plot point
419 get_s_plot(iplot, nplot, s);
420
421 // Get x position as Vector
422 interpolated_x(s, x);
423
424 // Get exact solution at this point
425 (*exact_soln_pt)(x, exact_soln);
426
427 // Output x,y,...,u_exact
428 for (unsigned i = 0; i < 2; i++)
429 {
430 outfile << x[i] << " ";
431 }
432 outfile << exact_soln[0] << std::endl;
433 }
434
435 // Write tecplot footer (e.g. FE connectivity lists)
436 write_tecplot_zone_footer(outfile, nplot);
437 }
438
439
440 //======================================================================
441 /// Validate against exact solution
442 ///
443 /// Solution is provided via function pointer.
444 /// Plot error at a given number of plot points.
445 ///
446 //======================================================================
448 std::ostream& outfile,
450 double& error,
451 double& norm)
452 {
453 // Initialise
454 error = 0.0;
455 norm = 0.0;
456
457 // Vector of local coordinates
458 Vector<double> s(2);
459
460 // Vector for coordintes
461 Vector<double> x(2);
462
463 // Find out how many nodes there are in the element
464 unsigned n_node = nnode();
465
466 Shape psi(n_node);
467
468 // Set the value of n_intpt
469 unsigned n_intpt = integral_pt()->nweight();
470
471 // Tecplot header info
472 outfile << "ZONE" << std::endl;
473
474 // Exact solution Vector (here a scalar)
475 Vector<double> exact_soln(1);
476
477 // Loop over the integration points
478 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
479 {
480 // Assign values of s
481 for (unsigned i = 0; i < 2; i++)
482 {
483 s[i] = integral_pt()->knot(ipt, i);
484 }
485
486 // Get the integral weight
487 double w = integral_pt()->weight(ipt);
488
489 // Get jacobian of mapping
490 double J = J_eulerian(s);
491
492 // Premultiply the weights and the Jacobian
493 double W = w * J;
494
495 // Get x position as Vector
496 interpolated_x(s, x);
497
498 // Get FE function value
500
501 // Get exact solution at this point
502 (*exact_soln_pt)(x, exact_soln);
503
504 // Output x,y,...,error
505 for (unsigned i = 0; i < 2; i++)
506 {
507 outfile << x[i] << " ";
508 }
509 outfile << exact_soln[0] << " " << exact_soln[0] - u_fe << std::endl;
510
511 // Add to error and norm
512 norm += exact_soln[0] * exact_soln[0] * x[0] * W;
513 error += (exact_soln[0] - u_fe) * (exact_soln[0] - u_fe) * x[0] * W;
514 }
515 }
516
517 //======================================================================
518 /// Calculate the integrated value of the unknown over the element
519 ///
520 //======================================================================
522 {
523 // Initialise
524 double sum = 0.0;
525
526 // Find out how many nodes there are in the element
527 const unsigned n_node = this->nnode();
528
529 // Find the index at which the concentration is stored
530 const unsigned u_nodal_index = this->u_index_cons_axisym_adv_diff();
531
532 // Allocate memory for the shape functions
533 Shape psi(n_node);
534
535 // Set the value of n_intpt
536 const unsigned n_intpt = integral_pt()->nweight();
537
538 // Loop over the integration points
539 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
540 {
541 // Get the integral weight
542 const double w = integral_pt()->weight(ipt);
543
544 // Get the shape functions
545 this->shape_at_knot(ipt, psi);
546
547 // Calculate the concentration
548 double interpolated_u = 0.0;
549 for (unsigned l = 0; l < n_node; l++)
550 {
551 interpolated_u += this->nodal_value(l, u_nodal_index) * psi(l);
552 }
553
554 // calculate the r coordinate
555 double r = 0.0;
556 for (unsigned l = 0; l < n_node; l++)
557 {
558 r += this->nodal_position(l, 0) * psi(l);
559 }
560
561 // Get jacobian of mapping
562 const double J = J_eulerian_at_knot(ipt);
563
564 // Add the values to the sum
565 sum += interpolated_u * r * w * J;
566 }
567
568 // return the sum
569 return sum;
570 }
571
572
573 //======================================================================
574 // Set the data for the number of Variables at each node
575 //======================================================================
576 template<unsigned NNODE_1D>
577 const unsigned
579
580 //====================================================================
581 // Force build of templates
582 //====================================================================
586
587} // 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
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
virtual double J_eulerian(const Vector< double > &s) const
Return the Jacobian of mapping from local to global coordinates at local position s.
Definition: elements.cc:4103
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
Definition: elements.h:2593
virtual std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction")
Definition: elements.h:3161
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3962
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
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as .
Definition: elements.h:1759
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &shifted_to_interior=false) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: elements.h:3148
virtual unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction")
Definition: elements.h:3186
double nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. If the node is hanging, the appropriate interpolation is ...
Definition: elements.h:2317
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
virtual double J_eulerian_at_knot(const unsigned &ipt) const
Return the Jacobian of the mapping from local to global coordinates at the ipt-th integration point.
Definition: elements.cc:4168
double raw_dnodal_position_dt(const unsigned &n, const unsigned &i) const
Return the i-th component of nodal velocity: dx/dt at local node n. Do not use the hanging node repre...
Definition: elements.h:2256
double raw_nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n but do NOT take hanging nodes into account.
Definition: elements.h:2576
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
virtual void write_tecplot_zone_footer(std::ostream &outfile, const unsigned &nplot) const
Add tecplot zone "footer" to output stream (when plotting nplot points in each "coordinate direction"...
Definition: elements.h:3174
virtual void shape_at_knot(const unsigned &ipt, Shape &psi) const
Return the geometric shape function at the ipt-th integration point.
Definition: elements.cc:3220
virtual unsigned self_test()
Self-test: Check inversion of element & do self-test for GeneralisedElement. Return 0 if OK.
Definition: elements.cc:4440
const double & pe_st() const
Peclet number multiplied by Strouhal number.
double interpolated_u_cons_axisym_adv_diff(const Vector< double > &s) const
Return FE representation of function value u(s) at local coordinate s.
double du_dt_cons_axisym_adv_diff(const unsigned &n) const
du/dt at local node n. Uses suitably interpolated value for hanging nodes.
virtual void get_conserved_wind_cons_axisym_adv_diff(const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, Vector< double > &wind) const
Get additional (conservative) wind at (Eulerian) position x and/or local coordinate s....
virtual double dshape_and_dtest_eulerian_at_knot_cons_axisym_adv_diff(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
Shape/test functions and derivs w.r.t. to global coords at integration point ipt; return Jacobian of ...
virtual void get_diff_cons_axisym_adv_diff(const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, DenseMatrix< double > &D) const
Get diffusivity tensor at (Eulerian) position x and/or local coordinate s. This function is virtual t...
static double Default_peclet_number
Static default value for the Peclet number.
virtual void get_wind_cons_axisym_adv_diff(const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, Vector< double > &wind) const
Get wind at (Eulerian) position x and/or local coordinate s. This function is virtual to allow overlo...
void output_fct(std::ostream &outfile, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output exact soln: r,z,u_exact at nplot^2 plot points.
double integrate_u()
Integrate the concentration over the element.
virtual void fill_in_generic_residual_contribution_cons_axisym_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.
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when time-derivatives are computed....
Function pointer to wind function Vector< double > & wind
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
virtual void get_source_cons_axisym_adv_diff(const unsigned &ipt, const Vector< double > &x, double &source) const
Get source term at (Eulerian) position x. This function is virtual to allow overloading in multi-phys...
void output(std::ostream &outfile)
Output with default number of plot points.
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
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.
//////////////////////////////////////////////////////////////////////// ////////////////////////////...
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...