womersley_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// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 
31 //Non-inline functions for Womersley elements
32 #include "womersley_elements.h"
33 
34 
35 namespace oomph
36 {
37 
38 
39 /// Default Womersley number
40 template<unsigned DIM>
41 double WomersleyEquations<DIM>::Default_ReSt_value = 0.0;
42 
43 
44 
45 //======================================================================
46 // Set the data for the number of Variables at each node
47 //======================================================================
48 template<unsigned DIM, unsigned NNODE_1D>
49 const unsigned QWomersleyElement<DIM,NNODE_1D>::Initial_Nvalue = 1;
50 
51 
52 //======================================================================
53 /// Compute element residual Vector and/or element Jacobian matrix
54 ///
55 /// flag=1: compute both
56 /// flag=0: compute only residual Vector
57 ///
58 /// Pure version without hanging nodes
59 //======================================================================
60 template <unsigned DIM>
63  DenseMatrix<double> &jacobian,
64  unsigned flag)
65 {
66  //Find out how many nodes there are
67  unsigned n_node = nnode();
68 
69  //Find the index at which the variable is stored
70  unsigned u_nodal_index = u_index_womersley();
71 
72  //Set up memory for the shape and test functions
73  Shape psi(n_node), test(n_node);
74  DShape dpsidx(n_node,DIM), dtestdx(n_node,DIM);
75 
76  //Set the value of n_intpt
77  unsigned n_intpt = integral_pt()->nweight();
78 
79  //Set the Vector to hold local coordinates
80  Vector<double> s(DIM);
81 
82  //Integers to hold the local equation and unknowns
83  int local_eqn=0, local_unknown=0;
84 
85  //Loop over the integration points
86  for(unsigned ipt=0;ipt<n_intpt;ipt++)
87  {
88 
89  //Assign values of s
90  for(unsigned i=0;i<DIM;i++) s[i] = integral_pt()->knot(ipt,i);
91 
92  //Get the integral weight
93  double w = integral_pt()->weight(ipt);
94 
95  //Call the derivatives of the shape and test functions
96  double J =
97  dshape_and_dtest_eulerian_at_knot_womersley(ipt,psi,dpsidx,test,dtestdx);
98 
99  //Premultiply the weights and the Jacobian
100  double W = w*J;
101 
102  //Allocate memory for local quantities and initialise to zero
103  double interpolated_u=0.0;
104  double dudt=0.0;
105  Vector<double> interpolated_x(DIM,0.0);
106  Vector<double> interpolated_dudx(DIM,0.0);
107 
108  //Calculate function value and derivatives:
109  // Loop over nodes
110  for(unsigned l=0;l<n_node;l++)
111  {
112  //Calculate the value at the nodes
113  double u_value = raw_nodal_value(l,u_nodal_index);
114  interpolated_u += u_value*psi(l);
115  dudt += du_dt_womersley(l)*psi(l);
116  // Loop over directions
117  for(unsigned j=0;j<DIM;j++)
118  {
119  interpolated_x[j] += raw_nodal_position(l,j)*psi(l);
120  interpolated_dudx[j] += u_value*dpsidx(l,j);
121  }
122  }
123 
124  //Get pressure gradient
125  //---------------------
126  double dpdz;
127 
128  //If no pressure gradient has been set, use zero
129  if(Pressure_gradient_data_pt==0)
130  {
131  dpdz = 0.0;
132  }
133  else
134  {
135  // Get gressure gradient
136  dpdz=Pressure_gradient_data_pt->value(0);
137  }
138 
139 
140  // Assemble residuals and Jacobian
141  //--------------------------------
142 
143  // Loop over the test functions
144  for(unsigned l=0;l<n_node;l++)
145  {
146  local_eqn = nodal_local_eqn(l,u_nodal_index);
147  /*IF it's not a boundary condition*/
148  if(local_eqn >= 0)
149  {
150 
151  // Add dpdz term and time derivative
152  residuals[local_eqn] += (re_st()*dudt+dpdz)*test(l)*W;
153 
154  // Laplace operator
155  for(unsigned k=0;k<DIM;k++)
156  {
157  residuals[local_eqn] +=
158  interpolated_dudx[k]*dtestdx(l,k)*W;
159  }
160 
161 
162  // Calculate the jacobian
163  //-----------------------
164  if(flag)
165  {
166  //Loop over the velocity shape functions again
167  for(unsigned l2=0;l2<n_node;l2++)
168  {
169  local_unknown = nodal_local_eqn(l2,u_nodal_index);
170  //If at a non-zero degree of freedom add in the entry
171  if(local_unknown >= 0)
172  {
173  // Mass matrix
174  jacobian(local_eqn,local_unknown)
175  += re_st()*test(l)*psi(l2)*
176  node_pt(l2)->time_stepper_pt()->weight(1,0)*W;
177 
178  // Laplace operator & mesh velocity bit
179  for(unsigned i=0;i<DIM;i++)
180  {
181  double tmp=dtestdx(l,i);
182  jacobian(local_eqn,local_unknown) += dpsidx(l2,i)*tmp*W;
183  }
184  }
185  }
186 
187  // Derivative w.r.t. pressure gradient data (if it's
188  // an unknown)
189  if ((Pressure_gradient_data_pt!=0)&&
190  (!Pressure_gradient_data_pt->is_pinned(0)))
191  {
192  local_unknown=external_local_eqn(0,0);
193  if (local_unknown>=0)
194  {
195  jacobian(local_eqn,local_unknown) += test(l)*W;
196  }
197 
198  // Derivatives of the final eqn (volume flux constraint
199  // w.r.t. to this unknown)
200  unsigned final_local_eqn=external_local_eqn(0,0);
201  unsigned local_unknown=local_eqn; // [from above just renaming
202  // // for clarity(!?)]
203  jacobian(final_local_eqn,local_unknown)+=psi(l)*W;
204  }
205 
206  }
207  }
208  }
209 
210 
211  } // End of loop over integration points
212 }
213 
214 
215 
216 
217 
218 //======================================================================
219 /// Compute volume flux through element
220 //======================================================================
221 template <unsigned DIM>
223 {
224  //Find out how many nodes there are
225  unsigned n_node = nnode();
226 
227  //Find the index at which the variable is stored
228  unsigned u_nodal_index = u_index_womersley();
229 
230  //Set up memory for the shape fcs
231  Shape psi(n_node);
232  DShape dpsidx(n_node,DIM);
233 
234  //Set the value of n_intpt
235  unsigned n_intpt = integral_pt()->nweight();
236 
237  //Set the Vector to hold local coordinates
238  Vector<double> s(DIM);
239 
240  // Initialise flux
241  double flux=0.0;
242 
243  //Loop over the integration points
244  for(unsigned ipt=0;ipt<n_intpt;ipt++)
245  {
246 
247  //Assign values of s
248  for(unsigned i=0;i<DIM;i++) s[i] = integral_pt()->knot(ipt,i);
249 
250  //Get the integral weight
251  double w = integral_pt()->weight(ipt);
252 
253  //Call the derivatives of the shape and test functions
254  double J = dshape_eulerian_at_knot(ipt,psi,dpsidx);
255 
256  //Premultiply the weights and the Jacobian
257  double W = w*J;
258 
259  //Allocate memory for local quantities and initialise to zero
260  double interpolated_u=0.0;
261 
262  //Calculate function value
263 
264  // Loop over nodes
265  for(unsigned l=0;l<n_node;l++)
266  {
267  // Calculate the value at the nodes (takes hanging node status
268  // into account
269  interpolated_u += nodal_value(l,u_nodal_index)*psi(l);
270  }
271 
272  // Add to flux
273  flux+=interpolated_u*W;
274 
275  } // End of loop over integration points
276 
277  return flux;
278 
279 }
280 
281 
282 
283 
284 
285 
286 //======================================================================
287 /// Self-test: Return 0 for OK
288 //======================================================================
289 template <unsigned DIM>
291 {
292 
293  bool passed=true;
294 
295  // Check lower-level stuff
296  if (FiniteElement::self_test()!=0)
297  {
298  passed=false;
299  }
300 
301  // Return verdict
302  if (passed)
303  {
304  return 0;
305  }
306  else
307  {
308  return 1;
309  }
310 
311 }
312 
313 
314 
315 
316 //======================================================================
317 /// Output function: x,y,z_out,0,0,u,0 to allow comparison
318 /// against full Navier Stokes at n_nplot x n_plot points (2D)
319 //======================================================================
320 template <unsigned DIM>
321 void WomersleyEquations<DIM>::output_3d(std::ostream &outfile,
322  const unsigned &nplot,
323  const double& z_out)
324 {
325  //Vector of local coordinates
326  Vector<double> s(DIM);
327 
328  // Tecplot header info
329  outfile << tecplot_zone_string(nplot);
330 
331  // Loop over plot points
332  unsigned num_plot_points=nplot_points(nplot);
333  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
334  {
335  // Get local coordinates of plot point
336  get_s_plot(iplot,nplot,s);
337 
338  for(unsigned i=0;i<DIM;i++)
339  {
340  outfile << interpolated_x(s,i) << " ";
341  }
342  outfile << z_out << " 0.0 0.0 " ;
343  outfile << interpolated_u_womersley(s);
344  outfile << " 0.0 " << std::endl;
345  }
346 
347  // Write tecplot footer (e.g. FE connectivity lists)
348  write_tecplot_zone_footer(outfile,nplot);
349 
350 }
351 
352 //======================================================================
353 /// Output function:
354 ///
355 /// x,y,u or x,y,z,u
356 ///
357 /// nplot points in each coordinate direction
358 //======================================================================
359 template <unsigned DIM>
360 void WomersleyEquations<DIM>::output(std::ostream &outfile,
361  const unsigned &nplot)
362 {
363  //Vector of local coordinates
364  Vector<double> s(DIM);
365 
366  // Tecplot header info
367  outfile << tecplot_zone_string(nplot);
368 
369  // Loop over plot points
370  unsigned num_plot_points=nplot_points(nplot);
371  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
372  {
373  // Get local coordinates of plot point
374  get_s_plot(iplot,nplot,s);
375 
376  for(unsigned i=0;i<DIM;i++)
377  {
378  outfile << interpolated_x(s,i) << " ";
379  }
380  outfile << interpolated_u_womersley(s) << 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 
390 //======================================================================
391 /// C-style output function:
392 ///
393 /// x,y,u or x,y,z,u
394 ///
395 /// nplot points in each coordinate direction
396 //======================================================================
397 template <unsigned DIM>
399  const unsigned &nplot)
400 {
401  //Vector of local coordinates
402  Vector<double> s(DIM);
403 
404  // Tecplot header info
405  fprintf(file_pt,"%s",tecplot_zone_string(nplot).c_str());
406 
407  // Loop over plot points
408  unsigned num_plot_points=nplot_points(nplot);
409  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
410  {
411 
412  // Get local coordinates of plot point
413  get_s_plot(iplot,nplot,s);
414 
415  for(unsigned i=0;i<DIM;i++)
416  {
417  fprintf(file_pt,"%g ",interpolated_x(s,i));
418 
419  }
420  fprintf(file_pt,"%g \n",interpolated_u_womersley(s));
421  }
422 
423  // Write tecplot footer (e.g. FE connectivity lists)
424  write_tecplot_zone_footer(file_pt,nplot);
425 
426 }
427 
428 
429 
430 //======================================================================
431  /// Output exact solution
432  ///
433  /// Solution is provided via function pointer.
434  /// Plot at a given number of plot points.
435  ///
436  /// x,y,u_exact or x,y,z,u_exact
437 //======================================================================
438 template <unsigned DIM>
439 void WomersleyEquations<DIM>::output_fct(std::ostream &outfile,
440  const unsigned &nplot,
442 {
443 
444  //Vector of local coordinates
445  Vector<double> s(DIM);
446 
447  // Vector for coordintes
448  Vector<double> x(DIM);
449 
450  // Tecplot header info
451  outfile << tecplot_zone_string(nplot);
452 
453  // Exact solution Vector (here a scalar)
454  Vector<double> exact_soln(1);
455 
456  // Loop over plot points
457  unsigned num_plot_points=nplot_points(nplot);
458  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
459  {
460 
461  // Get local coordinates of plot point
462  get_s_plot(iplot,nplot,s);
463 
464  // Get x position as Vector
465  interpolated_x(s,x);
466 
467  // Get exact solution at this point
468  (*exact_soln_pt)(x,exact_soln);
469 
470  //Output x,y,...,u_exact
471  for(unsigned i=0;i<DIM;i++)
472  {
473  outfile << x[i] << " ";
474  }
475  outfile << exact_soln[0] << std::endl;
476 
477  }
478 
479  // Write tecplot footer (e.g. FE connectivity lists)
480  write_tecplot_zone_footer(outfile,nplot);
481 
482 }
483 
484 
485 
486 //======================================================================
487 /// Output exact solution at time t
488 ///
489 /// Solution is provided via function pointer.
490 /// Plot at a given number of plot points.
491 ///
492 /// x,y,u_exact or x,y,z,u_exact
493 //======================================================================
494 template <unsigned DIM>
495 void WomersleyEquations<DIM>::output_fct(std::ostream &outfile,
496  const unsigned &nplot,
497  const double& time,
499 
500 {
501 
502  //Vector of local coordinates
503  Vector<double> s(DIM);
504 
505  // Vector for coordintes
506  Vector<double> x(DIM);
507 
508 
509  // Tecplot header info
510  outfile << tecplot_zone_string(nplot);
511 
512  // Exact solution Vector (here a scalar)
513  Vector<double> exact_soln(1);
514 
515  // Loop over plot points
516  unsigned num_plot_points=nplot_points(nplot);
517  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
518  {
519 
520  // Get local coordinates of plot point
521  get_s_plot(iplot,nplot,s);
522 
523  // Get x position as Vector
524  interpolated_x(s,x);
525 
526  // Get exact solution at this point
527  (*exact_soln_pt)(time,x,exact_soln);
528 
529  //Output x,y,...,u_exact
530  for(unsigned i=0;i<DIM;i++)
531  {
532  outfile << x[i] << " ";
533  }
534  outfile << exact_soln[0] << std::endl;
535 
536  }
537 
538  // Write tecplot footer (e.g. FE connectivity lists)
539  write_tecplot_zone_footer(outfile,nplot);
540 
541 }
542 
543 
544 
545 //======================================================================
546  /// Validate against exact solution
547  ///
548  /// Solution is provided via function pointer.
549  /// Plot error at a given number of plot points.
550  ///
551 //======================================================================
552 template <unsigned DIM>
554 compute_error(std::ostream &outfile,
556  double& error, double& norm)
557 {
558 
559  // Initialise
560  error=0.0;
561  norm=0.0;
562 
563  //Vector of local coordinates
564  Vector<double> s(DIM);
565 
566  // Vector for coordintes
567  Vector<double> x(DIM);
568 
569  //Find out how many nodes there are in the element
570  unsigned n_node = nnode();
571 
572  Shape psi(n_node);
573 
574  //Set the value of n_intpt
575  unsigned n_intpt = integral_pt()->nweight();
576 
577  // Tecplot header info
578  outfile << "ZONE" << std::endl;
579 
580  // Exact solution Vector (here a scalar)
581  Vector<double> exact_soln(1);
582 
583  //Loop over the integration points
584  for(unsigned ipt=0;ipt<n_intpt;ipt++)
585  {
586 
587  //Assign values of s
588  for(unsigned i=0;i<DIM;i++)
589  {
590  s[i] = integral_pt()->knot(ipt,i);
591  }
592 
593  //Get the integral weight
594  double w = integral_pt()->weight(ipt);
595 
596  // Get jacobian of mapping
597  double J=J_eulerian(s);
598 
599  //Premultiply the weights and the Jacobian
600  double W = w*J;
601 
602  // Get x position as Vector
603  interpolated_x(s,x);
604 
605  // Get FE function value
606  double u_fe = interpolated_u_womersley(s);
607 
608  // Get exact solution at this point
609  (*exact_soln_pt)(x,exact_soln);
610 
611  //Output x,y,...,error
612  for(unsigned i=0;i<DIM;i++)
613  {
614  outfile << x[i] << " ";
615  }
616  outfile << exact_soln[0] << " " << exact_soln[0]-u_fe << std::endl;
617 
618  // Add to error and norm
619  norm+=exact_soln[0]*exact_soln[0]*W;
620  error+=(exact_soln[0]-u_fe)*(exact_soln[0]-u_fe)*W;
621  }
622 }
623 
624 
625 
626 
627 //======================================================================
628 /// Validate against exact solution at time t.
629 ///
630 /// Solution is provided via function pointer.
631 /// Plot error at a given number of plot points.
632 ///
633 //======================================================================
634 template<unsigned DIM>
635 void WomersleyEquations<DIM>::compute_error(std::ostream &outfile,
637  const double& time, double& error, double& norm)
638 
639 {
640 
641  // Initialise
642  error=0.0;
643  norm=0.0;
644 
645  //Vector of local coordinates
646  Vector<double> s(DIM);
647 
648  // Vector for coordintes
649  Vector<double> x(DIM);
650 
651 
652  //Find out how many nodes there are in the element
653  unsigned n_node = nnode();
654 
655  Shape psi(n_node);
656 
657  //Set the value of n_intpt
658  unsigned n_intpt = integral_pt()->nweight();
659 
660  // Tecplot
661  outfile << "ZONE" << std::endl;
662 
663  // Exact solution Vector (here a scalar)
664  Vector<double> exact_soln(1);
665 
666  //Loop over the integration points
667  for(unsigned ipt=0;ipt<n_intpt;ipt++)
668  {
669 
670  //Assign values of s
671  for(unsigned i=0;i<DIM;i++)
672  {
673  s[i] = integral_pt()->knot(ipt,i);
674  }
675 
676  //Get the integral weight
677  double w = integral_pt()->weight(ipt);
678 
679  // Get jacobian of mapping
680  double J=J_eulerian(s);
681 
682  //Premultiply the weights and the Jacobian
683  double W = w*J;
684 
685  // Get x position as Vector
686  interpolated_x(s,x);
687 
688  // Get FE function value
689  double u_fe = interpolated_u_womersley(s);
690 
691  // Get exact solution at this point
692  (*exact_soln_pt)(time,x,exact_soln);
693 
694  //Output x,y,...,error
695  for(unsigned i=0;i<DIM;i++)
696  {
697  outfile << x[i] << " ";
698  }
699  outfile << exact_soln[0] << " " << exact_soln[0]-u_fe << std::endl;
700 
701  // Add to error and norm
702  norm+=exact_soln[0]*exact_soln[0]*W;
703  error+=(exact_soln[0]-u_fe)*(exact_soln[0]-u_fe)*W;
704 
705  }
706 }
707 
708 
709 
710 
711 //====================================================================
712 // Force build of templates
713 //====================================================================
714 
715 template class WomersleyEquations<1>;
716 
717 template class QWomersleyElement<1,2>;
718 template class QWomersleyElement<1,3>;
719 template class QWomersleyElement<1,4>;
720 
721 template class WomersleyEquations<2>;
722 
723 template class QWomersleyElement<2,2>;
724 template class QWomersleyElement<2,3>;
725 template class QWomersleyElement<2,4>;
726 
727 //template class QWomersleyElement<3,2>;
728 //template class QWomersleyElement<3,3>;
729 //template class QWomersleyElement<3,4>;
730 
731 }
732 
unsigned self_test()
Self-test: Return 0 for OK.
virtual void fill_in_generic_residual_contribution_womersley(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
Compute element residual Vector only (if flag=and/or element Jacobian matrix.
cstr elem_len * i
Definition: cfortran.h:607
virtual unsigned self_test()
Self-test: Check inversion of element & do self-test for GeneralisedElement. Return 0 if OK...
Definition: elements.cc:4247
double get_volume_flux()
Compute total volume flux through element.
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
static char t char * s
Definition: cfortran.h:572
void output_3d(std::ostream &outfile, const unsigned &n_plot, const double &z_out)
Output function: x,y,z_out,0,0,u,0 to allow comparison against full Navier Stokes at n_nplot x n_plot...
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
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.
void output(std::ostream &outfile)
Output with default number of plot points.
void(* UnsteadyExactSolutionFctPt)(const double &, const Vector< double > &, Vector< double > &)
Function pointer for function that computes Vector-valued time-dependent function as ...
Definition: elements.h:1726