linear_wave_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 #include <cmath>
31 
32 //Non-inline functions for LinearWave elements
33 #include "linear_wave_elements.h"
34 
35 
36 namespace oomph
37 {
38 
39 
40 ///2D LinearWave elements
41 
42 
43 //======================================================================
44 // Set the data for the number of Variables at each node
45 //======================================================================
46 template<>
47 const unsigned QLinearWaveElement<1,4>::Initial_Nvalue[4]={1,1,1,1};
48 template<>
49 const unsigned QLinearWaveElement<1,3>::Initial_Nvalue[3]={1,1,1};
50 template<>
52 
53 template<>
55 ={1,1,1,1,1,1,1,1,
56  1,1,1,1,1,1,1,1};
57 template<>
59 ={1,1,1,1,1,1,1,1,1};
60 template<>
62 ={1,1,1,1};
63 
64 template<>
66 ={1,1,1,1,1,1,1,1};
67 template<>
69 ={1,1,1,1,1,1,1,1,1,
70  1,1,1,1,1,1,1,1,1,
71  1,1,1,1,1,1,1,1,1};
72 template<>
74 ={1,1,1,1,1,1,1,1,
75  1,1,1,1,1,1,1,1,
76  1,1,1,1,1,1,1,1,
77  1,1,1,1,1,1,1,1,
78  1,1,1,1,1,1,1,1,
79  1,1,1,1,1,1,1,1,
80  1,1,1,1,1,1,1,1,
81  1,1,1,1,1,1,1,1};
82 
83 
84 
85 //======================================================================
86 /// Compute element residual Vector and/or element Jacobian matrix
87 ///
88 /// flag=1: compute both
89 /// flag=0: compute only residual Vector
90 ///
91 /// Pure version without hanging nodes
92 //======================================================================
93 template <unsigned DIM>
96  DenseMatrix<double> &jacobian,
97  unsigned flag)
98 {
99  //Find out how many nodes there are
100  unsigned n_node = nnode();
101 
102  // Get continuous time from timestepper of first node
103  double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
104 
105  // Find the index at which the linear wave unknown is stored
106  unsigned u_nodal_index = u_index_lin_wave();
107 
108  //Set up memory for the shape and test functions
109  Shape psi(n_node), test(n_node);
110  DShape dpsidx(n_node,DIM), dtestdx(n_node,DIM);
111 
112  //Set the value of n_intpt
113  unsigned n_intpt = integral_pt()->nweight();
114 
115  //Set the Vector to hold local coordinates
116  Vector<double> s(DIM);
117 
118  //Integers to hold the local equation and unknowns
119  int local_eqn=0, local_unknown=0;
120 
121  //Loop over the integration points
122  for(unsigned ipt=0;ipt<n_intpt;ipt++)
123  {
124 
125  //Assign values of s
126  for(unsigned i=0;i<DIM;i++) s[i] = integral_pt()->knot(ipt,i);
127 
128  //Get the integral weight
129  double w = integral_pt()->weight(ipt);
130 
131  //Call the derivatives of the shape and test functions
132  double J =
133  dshape_and_dtest_eulerian_at_knot_lin_wave(ipt,psi,dpsidx,test,dtestdx);
134 
135  //Premultiply the weights and the Jacobian
136  double W = w*J;
137 
138  //Allocate memory for local quantities and initialise to zero
139  double interpolated_u=0.0;
140  double ddudt=0.0;
141  Vector<double> interpolated_x(DIM,0.0);
142  Vector<double> interpolated_dudx(DIM,0.0);
143 
144  //Calculate function value and derivatives:
145  // Loop over nodes
146  for(unsigned l=0;l<n_node;l++)
147  {
148  //Get the nodal value
149  double u_value = raw_nodal_value(l,u_nodal_index);
150  interpolated_u += u_value*psi(l);
151  ddudt += d2u_dt2_lin_wave(l)*psi(l);
152  // Loop over directions
153  for(unsigned j=0;j<DIM;j++)
154  {
155  interpolated_x[j] += raw_nodal_position(l,j)*psi(l);
156  interpolated_dudx[j] += u_value*dpsidx(l,j);
157  }
158  }
159 
160 
161  //Get source function
162  //-------------------
163  double source;
164  get_source_lin_wave(time,ipt,interpolated_x,source);
165 
166  // Assemble residuals and Jacobian
167  //--------------------------------
168 
169  // Loop over the test functions
170  for(unsigned l=0;l<n_node;l++)
171  {
172  local_eqn = nodal_local_eqn(l,u_nodal_index);
173  /*IF it's not a boundary condition*/
174  if(local_eqn >= 0)
175  {
176 
177  // Add body force/source term and time derivative
178  residuals[local_eqn] += (ddudt+source)*test(l)*W;
179 
180  // Laplace operator
181  for(unsigned k=0;k<DIM;k++)
182  {
183  residuals[local_eqn] +=
184  interpolated_dudx[k]*dtestdx(l,k)*W;
185  }
186 
187 
188  // Calculate the jacobian
189  //-----------------------
190  if(flag)
191  {
192  //Loop over the velocity shape functions again
193  for(unsigned l2=0;l2<n_node;l2++)
194  {
195  local_unknown = nodal_local_eqn(l2,u_nodal_index);
196  //If at a non-zero degree of freedom add in the entry
197  if(local_unknown >= 0)
198  {
199  // Mass matrix
200  jacobian(local_eqn,local_unknown)
201  += test(l)*psi(l2)*node_pt(l2)->time_stepper_pt()->weight(2,0)*W;
202 
203  // Laplace operator
204  for(unsigned i=0;i<DIM;i++)
205  {
206  jacobian(local_eqn,local_unknown)
207  += dpsidx(l2,i)*dtestdx(l,i)*W;
208  }
209  }
210  }
211  }
212  }
213  }
214 
215 
216  } // End of loop over integration points
217 }
218 
219 
220 
221 
222 //======================================================================
223 /// Self-test: Return 0 for OK
224 //======================================================================
225 template <unsigned DIM>
227 {
228 
229  bool passed=true;
230 
231  // Check lower-level stuff
232  if (FiniteElement::self_test()!=0)
233  {
234  passed=false;
235  }
236 
237  // Return verdict
238  if (passed)
239  {
240  return 0;
241  }
242  else
243  {
244  return 1;
245  }
246 
247 }
248 
249 
250 
251 //======================================================================
252 /// Output function:
253 ///
254 /// x,y,u or x,y,z,u
255 ///
256 /// nplot points in each coordinate direction
257 //======================================================================
258 template <unsigned DIM>
259 void LinearWaveEquations<DIM>::output(std::ostream &outfile,
260  const unsigned &nplot)
261 {
262  //Vector of local coordinates
263  Vector<double> s(DIM);
264 
265  // Tecplot header info
266  outfile << tecplot_zone_string(nplot);
267 
268  // Loop over plot points
269  unsigned num_plot_points=nplot_points(nplot);
270  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
271  {
272  // Get local coordinates of plot point
273  get_s_plot(iplot,nplot,s);
274  for(unsigned i=0;i<DIM;i++)
275  {
276  outfile << interpolated_x(s,i) << " ";
277  }
278  outfile << interpolated_u_lin_wave(s) << " ";
279  outfile << interpolated_du_dt_lin_wave(s) << " ";
280  outfile << interpolated_d2u_dt2_lin_wave(s) << std::endl;
281  }
282 
283  // Write tecplot footer (e.g. FE connectivity lists)
284  write_tecplot_zone_footer(outfile,nplot);
285 
286 }
287 
288 //======================================================================
289 /// C-style output function:
290 ///
291 /// x,y,u or x,y,z,u
292 ///
293 /// nplot points in each coordinate direction
294 //======================================================================
295 template <unsigned DIM>
297  const unsigned &nplot)
298 {
299  //Vector of local coordinates
300  Vector<double> s(DIM);
301 
302  // Tecplot header info
303  fprintf(file_pt,"%s",tecplot_zone_string(nplot).c_str());
304 
305  // Loop over plot points
306  unsigned num_plot_points=nplot_points(nplot);
307  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
308  {
309  // Get local coordinates of plot point
310  get_s_plot(iplot,nplot,s);
311 
312  for(unsigned i=0;i<DIM;i++)
313  {
314  //outfile << interpolated_x(s,i) << " ";
315  fprintf(file_pt,"%g ",interpolated_x(s,i));
316  }
317  //outfile << interpolated_u(s) << std::endl;
318  fprintf(file_pt,"%g \n",interpolated_u_lin_wave(s));
319  }
320 
321  // Write tecplot footer (e.g. FE connectivity lists)
322  write_tecplot_zone_footer(file_pt,nplot);
323 
324 }
325 
326 
327 
328 //======================================================================
329 /// Output exact solution
330 ///
331 /// Solution is provided via function pointer.
332 /// Plot at a given number of plot points.
333 ///
334 /// x,y,u_exact or x,y,z,u_exact
335 //======================================================================
336 template <unsigned DIM>
338  std::ostream &outfile, const unsigned &nplot,
340 {
341  //Vector of local coordinates
342  Vector<double> s(DIM);
343 
344  // Vector for coordintes
345  Vector<double> x(DIM);
346 
347  // Tecplot header info
348  outfile << tecplot_zone_string(nplot);
349 
350  // Exact solution Vector (value and first and second time derivs)
351  Vector<double> exact_soln(3);
352 
353  // Loop over plot points
354  unsigned num_plot_points=nplot_points(nplot);
355  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
356  {
357  // Get local coordinates of plot point
358  get_s_plot(iplot,nplot,s);
359 
360  // Get x position as Vector
361  interpolated_x(s,x);
362 
363  // Get exact solution at this point
364  (*exact_soln_pt)(x,exact_soln);
365 
366  //Output x,y,...,u_exact
367  for(unsigned i=0;i<DIM;i++)
368  {
369  outfile << x[i] << " ";
370  }
371  outfile << exact_soln[0] << std::endl;
372 
373  }
374 
375  // Write tecplot footer (e.g. FE connectivity lists)
376  write_tecplot_zone_footer(outfile,nplot);
377 }
378 
379 
380 
381 //======================================================================
382 /// Output exact solution at time t
383 ///
384 /// Solution is provided via function pointer.
385 /// Plot at a given number of plot points.
386 ///
387 /// x,y,u_exact or x,y,z,u_exact
388 //======================================================================
389 template <unsigned DIM>
391 std::ostream &outfile, const unsigned &nplot, const double& time,
393 {
394 
395  //Vector of local coordinates
396  Vector<double> s(DIM);
397 
398  // Vector for coordintes
399  Vector<double> x(DIM);
400 
401  // Tecplot header info
402  outfile << tecplot_zone_string(nplot);
403 
404  // Exact solution Vector u, dudt, d2udt2
405  Vector<double> exact_soln(3);
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  // Get x position as Vector
416  interpolated_x(s,x);
417 
418  // Get exact solution at this point
419  (*exact_soln_pt)(time,x,exact_soln);
420 
421  //Output x,y,...,u_exact
422  for(unsigned i=0;i<DIM;i++)
423  {
424  outfile << x[i] << " ";
425  }
426  outfile << exact_soln[0] << " "
427  << exact_soln[1] << " "
428  << exact_soln[2] << " "
429  << std::endl;
430  }
431 
432  // Write tecplot footer (e.g. FE connectivity lists)
433  write_tecplot_zone_footer(outfile,nplot);
434 }
435 
436 
437 
438 
439 //======================================================================
440 /// Validate against exact solution
441 ///
442 /// Solution is provided via function pointer.
443 /// Plot error at a given number of plot points.
444 ///
445 //======================================================================
446 template <unsigned DIM>
448  std::ostream &outfile,
450  double& error, double& norm)
451 {
452 
453  // Initialise
454  error=0.0;
455  norm=0.0;
456 
457  //Vector of local coordinates
458  Vector<double> s(DIM);
459 
460  // Vector for coordintes
461  Vector<double> x(DIM);
462 
463 
464  //Find out how many nodes there are in the element
465  unsigned n_node = nnode();
466 
467  Shape psi(n_node);
468 
469  //Set the value of n_intpt
470  unsigned n_intpt = integral_pt()->nweight();
471 
472  // Tecplot header info
473  outfile << "ZONE" << std::endl;
474 
475  // Exact solution Vector (value, and first and second time derivs)
476  Vector<double> exact_soln(3);
477 
478  //Loop over the integration points
479  for(unsigned ipt=0;ipt<n_intpt;ipt++)
480  {
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_lin_wave(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 
521 
522 //======================================================================
523 /// Validate against exact solution at time t.
524 ///
525 /// Solution is provided via function pointer.
526 /// Plot error at a given number of plot points.
527 ///
528 //======================================================================
529 template <unsigned DIM>
531  std::ostream &outfile,
533  const double& time, double& error, double& norm)
534 
535 {
536 
537  // Initialise
538  error=0.0;
539  norm=0.0;
540 
541  //Vector of local coordinates
542  Vector<double> s(DIM);
543 
544  // Vector for coordintes
545  Vector<double> x(DIM);
546 
547  //Find out how many nodes there are in the element
548  unsigned n_node = nnode();
549 
550  Shape psi(n_node);
551 
552  //Set the value of n_intpt
553  unsigned n_intpt = integral_pt()->nweight();
554 
555  // Tecplot header info
556  outfile << "ZONE" << std::endl;
557 
558  // Exact solution Vector (value, and first and second time derivs)
559  Vector<double> exact_soln(3);
560 
561  //Loop over the integration points
562  for(unsigned ipt=0;ipt<n_intpt;ipt++)
563  {
564 
565  //Assign values of s
566  for(unsigned i=0;i<DIM;i++)
567  {
568  s[i] = integral_pt()->knot(ipt,i);
569  }
570 
571  //Get the integral weight
572  double w = integral_pt()->weight(ipt);
573 
574  // Get jacobian of mapping
575  double J=J_eulerian(s);
576 
577  //Premultiply the weights and the Jacobian
578  double W = w*J;
579 
580  // Get x position as Vector
581  interpolated_x(s,x);
582 
583  // Get FE function value
584  double u_fe=interpolated_u_lin_wave(s);
585 
586  // Get exact solution at this point
587  (*exact_soln_pt)(time,x,exact_soln);
588 
589  //Output x,y,...,error
590  for(unsigned i=0;i<DIM;i++)
591  {
592  outfile << x[i] << " ";
593  }
594  outfile << exact_soln[0] << " " << exact_soln[0]-u_fe << std::endl;
595 
596  // Add to error and norm
597  norm+=exact_soln[0]*exact_soln[0]*W;
598  error+=(exact_soln[0]-u_fe)*(exact_soln[0]-u_fe)*W;
599 
600  }
601 }
602 
603 
604 
605 
606 //====================================================================
607 // Force build of templates
608 //====================================================================
609 template class QLinearWaveElement<1,2>;
610 template class QLinearWaveElement<1,3>;
611 template class QLinearWaveElement<1,4>;
612 
613 template class QLinearWaveElement<2,2>;
614 template class QLinearWaveElement<2,3>;
615 template class QLinearWaveElement<2,4>;
616 
617 template class QLinearWaveElement<3,2>;
618 template class QLinearWaveElement<3,3>;
619 template class QLinearWaveElement<3,4>;
620 
621 
622 }
virtual void fill_in_generic_residual_contribution_lin_wave(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
static char t char * s
Definition: cfortran.h:572
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(std::ostream &outfile)
Output with default number of plot points.
unsigned self_test()
Self-test: Return 0 for OK.
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
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
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.