harmonic.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 //Driver to solve the harmonic equation with homogeneous Dirichlet boundary
31 //conditions.
32 
33 // Generic oomph-lib routines
34 #include "generic.h"
35 
36 // Include the mesh
37 #include "meshes/one_d_mesh.h"
38 
39 using namespace std;
40 
41 using namespace oomph;
42 
43 
44 //===================================================================
45 /// Function-type-object to perform comparison of complex data types
46 /// Needed to sort the complex eigenvalues into order based on the
47 /// size of the real part.
48 //==================================================================
49 template <class T>
50 class ComplexLess
51 {
52 public:
53 
54  /// Comparison. Are the values identical or not?
55  bool operator()(const complex<T> &x, const complex<T> &y) const
56  {
57  return x.real() < y.real();
58  }
59 };
60 
61 
62 //=================================================================
63 /// A class for all elements that solve the simple one-dimensional
64 /// eigenvalue problem
65 /// \f[
66 /// \frac{\partial^2 u}{\partial x_i^2} + \lambda u = 0
67 /// \f]
68 /// These elements are very closely related to the Poisson
69 /// elements and could inherit from them. They are here developed
70 /// from scratch for pedagogical purposes.
71 /// This class contains the generic maths. Shape functions, geometric
72 /// mapping etc. must get implemented in derived class.
73 //================================================================
74 class HarmonicEquations : public virtual FiniteElement
75 {
76 
77 public:
78  /// Empty Constructor
80 
81  /// \short Access function: Eigenfunction value at local node n
82  /// Note that solving the eigenproblem does not assign values
83  /// to this storage space. It is used for output purposes only.
84  virtual inline double u(const unsigned& n) const
85  {return nodal_value(n,0);}
86 
87  /// Output the eigenfunction with default number of plot points
88  void output(ostream &outfile)
89  {
90  unsigned nplot=5;
91  output(outfile,nplot);
92  }
93 
94  /// \short Output FE representation of soln: x,y,u or x,y,z,u at
95  /// Nplot plot points
96  void output(ostream &outfile, const unsigned &nplot)
97  {
98  //Vector of local coordinates
99  Vector<double> s(1);
100 
101  // Tecplot header info
102  outfile << tecplot_zone_string(nplot);
103 
104  // Loop over plot points
105  unsigned num_plot_points=nplot_points(nplot);
106  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
107  {
108  // Get local coordinates of plot point
109  get_s_plot(iplot,nplot,s);
110  //Output the coordinate and the eigenfunction
111  outfile << interpolated_x(s,0) << " " << interpolated_u(s) << std::endl;
112  }
113  // Write tecplot footer (e.g. FE connectivity lists)
114  write_tecplot_zone_footer(outfile,nplot);
115  }
116 
117  /// \short Assemble the contributions to the jacobian and mass
118  /// matrices
120  Vector<double> &residuals,
121  DenseMatrix<double> &jacobian, DenseMatrix<double> &mass_matrix)
122  {
123  //Find out how many nodes there are
124  unsigned n_node = nnode();
125 
126  //Set up memory for the shape functions and their derivatives
127  Shape psi(n_node);
128  DShape dpsidx(n_node,1);
129 
130  //Set the number of integration points
131  unsigned n_intpt = integral_pt()->nweight();
132 
133  //Integers to store the local equation and unknown numbers
134  int local_eqn=0, local_unknown=0;
135 
136  //Loop over the integration points
137  for(unsigned ipt=0;ipt<n_intpt;ipt++)
138  {
139  //Get the integral weight
140  double w = integral_pt()->weight(ipt);
141 
142  //Call the derivatives of the shape and test functions
143  double J = dshape_eulerian_at_knot(ipt,psi,dpsidx);
144 
145  //Premultiply the weights and the Jacobian
146  double W = w*J;
147 
148  //Assemble the contributions to the mass matrix
149  //Loop over the test functions
150  for(unsigned l=0;l<n_node;l++)
151  {
152  //Get the local equation number
153  local_eqn = u_local_eqn(l);
154  /*IF it's not a boundary condition*/
155  if(local_eqn >= 0)
156  {
157  //Loop over the shape functions
158  for(unsigned l2=0;l2<n_node;l2++)
159  {
160  local_unknown = u_local_eqn(l2);
161  //If at a non-zero degree of freedom add in the entry
162  if(local_unknown >= 0)
163  {
164  jacobian(local_eqn,local_unknown) += dpsidx(l,0)*dpsidx(l2,0)*W;
165  mass_matrix(local_eqn, local_unknown) += psi(l)*psi(l2)*W;
166  }
167  }
168  }
169  }
170  }
171  } //end_of_fill_in_contribution_to_jacobian_and_mass_matrix
172 
173  /// Return FE representation of function value u(s) at local coordinate s
174  inline double interpolated_u(const Vector<double> &s) const
175  {
176  unsigned n_node = nnode();
177 
178  //Local shape functions
179  Shape psi(n_node);
180 
181  //Find values of basis function
182  this->shape(s,psi);
183 
184  //Initialise value of u
185  double interpolated_u = 0.0;
186 
187  //Loop over the local nodes and sum
188  for(unsigned l=0;l<n_node;l++) {interpolated_u+=u(l)*psi[l];}
189 
190  //Return the interpolated value of the eigenfunction
191  return(interpolated_u);
192  }
193 
194 protected:
195 
196  /// \short Shape/test functions and derivs w.r.t. to global coords at
197  /// local coord. s; return Jacobian of mapping
198  virtual double dshape_eulerian(const Vector<double> &s,
199  Shape &psi,
200  DShape &dpsidx) const=0;
201 
202  /// \short Shape/test functions and derivs w.r.t. to global coords at
203  /// integration point ipt; return Jacobian of mapping
204  virtual double dshape_eulerian_at_knot(const unsigned &ipt,
205  Shape &psi,
206  DShape &dpsidx) const=0;
207 
208  /// \short Access function that returns the local equation number
209  /// of the unknown in the problem. Default is to assume that it is the
210  /// first (only) value stored at the node.
211  virtual inline int u_local_eqn(const unsigned &n)
212  {return nodal_local_eqn(n,0);}
213 
214  private:
215 
216 };
217 
218 
219 
220 //======================================================================
221 /// QHarmonicElement<NNODE_1D> elements are 1D Elements with
222 /// NNODE_1D nodal points that are used to solve the Harmonic eigenvalue
223 /// Problem described by HarmonicEquations.
224 //======================================================================
225 template <unsigned NNODE_1D>
226 class QHarmonicElement : public virtual QElement<1,NNODE_1D>,
227  public HarmonicEquations
228 {
229 
230  public:
231 
232  ///\short Constructor: Call constructors for QElement and
233  /// Poisson equations
234  QHarmonicElement() : QElement<1,NNODE_1D>(), HarmonicEquations() {}
235 
236  /// \short Required # of `values' (pinned or dofs)
237  /// at node n
238  inline unsigned required_nvalue(const unsigned &n) const {return 1;}
239 
240  /// \short Output function overloaded from HarmonicEquations
241  void output(ostream &outfile)
242  {HarmonicEquations::output(outfile);}
243 
244  /// \short Output function overloaded from HarmonicEquations
245  void output(ostream &outfile, const unsigned &Nplot)
246  {HarmonicEquations::output(outfile,Nplot);}
247 
248 
249 protected:
250 
251 /// Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
252  inline double dshape_eulerian(const Vector<double> &s,
253  Shape &psi,
254  DShape &dpsidx) const
255  {return QElement<1,NNODE_1D>::dshape_eulerian(s,psi,dpsidx);}
256 
257 
258  /// \short Shape, test functions & derivs. w.r.t. to global coords. at
259  /// integration point ipt. Return Jacobian.
260  inline double dshape_eulerian_at_knot(const unsigned& ipt,
261  Shape &psi,
262  DShape &dpsidx) const
263  {return QElement<1,NNODE_1D>::dshape_eulerian_at_knot(ipt,psi,dpsidx);}
264 
265 }; //end_of_QHarmonic_class_definition
266 
267 
268 //==start_of_problem_class============================================
269 /// 1D Harmonic problem in unit interval.
270 //====================================================================
271 template<class ELEMENT,class EIGEN_SOLVER>
272 class HarmonicProblem : public Problem
273 {
274 public:
275 
276  /// Constructor: Pass number of elements and pointer to source function
277  HarmonicProblem(const unsigned& n_element);
278 
279  /// Destructor (empty)
280  ~HarmonicProblem(){delete this->mesh_pt(); delete this->eigen_solver_pt();}
281 
282  /// Solve the problem
283  void solve(const unsigned &label);
284 
285  /// \short Doc the solution, pass the number of the case considered,
286  /// so that output files can be distinguished.
287  void doc_solution(const unsigned& label);
288 
289 }; // end of problem class
290 
291 
292 
293 //=====start_of_constructor===============================================
294 /// \short Constructor for 1D Harmonic problem in unit interval.
295 /// Discretise the 1D domain with n_element elements of type ELEMENT.
296 /// Specify function pointer to source function.
297 //========================================================================
298 template<class ELEMENT,class EIGEN_SOLVER>
300  const unsigned& n_element)
301 {
302  //Create the eigen solver
303  this->eigen_solver_pt() = new EIGEN_SOLVER;
304 
305  //Get the positive eigenvalues, shift is zero by default
306  static_cast<EIGEN_SOLVER*>(eigen_solver_pt())
307  ->get_eigenvalues_right_of_shift();
308 
309  //Set domain length
310  double L=1.0;
311 
312  // Build mesh and store pointer in Problem
313  Problem::mesh_pt() = new OneDMesh<ELEMENT>(n_element,L);
314 
315  // Set the boundary conditions for this problem: By default, all nodal
316  // values are free -- we only need to pin the ones that have
317  // Dirichlet conditions.
318 
319  // Pin the single nodal value at the single node on mesh
320  // boundary 0 (= the left domain boundary at x=0)
321  mesh_pt()->boundary_node_pt(0,0)->pin(0);
322 
323  // Pin the single nodal value at the single node on mesh
324  // boundary 1 (= the right domain boundary at x=1)
325  mesh_pt()->boundary_node_pt(1,0)->pin(0);
326 
327  // Setup equation numbering scheme
328  assign_eqn_numbers();
329 
330 } // end of constructor
331 
332 
333 
334 //===start_of_doc=========================================================
335 /// Doc the solution in tecplot format. Label files with label.
336 //========================================================================
337 template<class ELEMENT,class EIGEN_SOLVER>
339 {
340 
341  ofstream some_file;
342  char filename[100];
343 
344  // Number of plot points
345  unsigned npts;
346  npts=5;
347 
348  // Output solution with specified number of plot points per element
349  sprintf(filename,"soln%i.dat",label);
350  some_file.open(filename);
351  mesh_pt()->output(some_file,npts);
352  some_file.close();
353 
354 } // end of doc
355 
356 //=======================start_of_solve==============================
357 /// Solve the eigenproblem
358 //===================================================================
359 template<class ELEMENT,class EIGEN_SOLVER>
361 solve(const unsigned& label)
362 {
363  //Set external storage for the eigenvalues
364  Vector<complex<double> > eigenvalues;
365  //Set external storage for the eigenvectors
366  Vector<DoubleVector> eigenvectors;
367  //Desired number eigenvalues
368  unsigned n_eval=4;
369 
370  //Solve the eigenproblem
371  this->solve_eigenproblem(n_eval,eigenvalues,eigenvectors);
372 
373  //We now need to sort the output based on the size of the real part
374  //of the eigenvalues.
375  //This is because the solver does not necessarily sort the eigenvalues
376  Vector<complex<double> > sorted_eigenvalues = eigenvalues;
377  sort(sorted_eigenvalues.begin(),sorted_eigenvalues.end(),
379 
380  //Read out the second smallest eigenvalue
381  complex<double> temp_evalue = sorted_eigenvalues[1];
382  unsigned second_smallest_index=0;
383  //Loop over the unsorted eigenvalues and find the entry that corresponds
384  //to our second smallest eigenvalue.
385  for(unsigned i=0;i<eigenvalues.size();i++)
386  {
387  //Note that equality tests for doubles are bad, but it was just
388  //sorted data, so should be fine
389  if(eigenvalues[i] == temp_evalue) {second_smallest_index=i; break;}
390  }
391 
392  //Normalise the eigenvector
393  {
394  //Get the dimension of the eigenvector
395  unsigned dim = eigenvectors[second_smallest_index].nrow();
396  double length=0.0;
397  //Loop over all the entries
398  for(unsigned i=0;i<dim;i++)
399  {
400  //Add the contribution to the length
401  length += std::pow(eigenvectors[second_smallest_index][i],2.0);
402  }
403  //Now take the magnitude
404  length = sqrt(length);
405  //Fix the sign
406  if(eigenvectors[second_smallest_index][0] < 0) {length *= -1.0;}
407  //Finally normalise
408  for(unsigned i=0;i<dim;i++)
409  {
410  eigenvectors[second_smallest_index][i] /= length;
411  }
412  }
413 
414  //Now assign the second eigenvector to the dofs of the problem
415  this->assign_eigenvector_to_dofs(eigenvectors[second_smallest_index]);
416  //Output solution for this case (label output files with "1")
417  this->doc_solution(label);
418 
419  char filename[100];
420  sprintf(filename,"eigenvalues%i.dat",label);
421 
422  //Open an output file for the sorted eigenvalues
423  ofstream evalues(filename);
424  for(unsigned i=0;i<n_eval;i++)
425  {
426  //Print to screen
427  cout << sorted_eigenvalues[i].real() << " "
428  << sorted_eigenvalues[i].imag() << std::endl;
429  //Send to file
430  evalues << sorted_eigenvalues[i].real() << " "
431  << sorted_eigenvalues[i].imag() << std::endl;
432  }
433 
434  evalues.close();
435 } //end_of_solve
436 
437 
438 ////////////////////////////////////////////////////////////////////////
439 ////////////////////////////////////////////////////////////////////////
440 ////////////////////////////////////////////////////////////////////////
441 
442 
443 //======start_of_main==================================================
444 /// Driver for 1D Poisson problem
445 //=====================================================================
446 int main(int argc, char **argv)
447 {
448 //Want to test Trilinos if we have it, so we must initialise MPI
449 //if we have compiled with it
450 #ifdef OOMPH_HAS_MPI
451  MPI_Helpers::init(argc,argv);
452 #endif
453 
454  // Set up the problem:
455  unsigned n_element=100; //Number of elements
456 
457  clock_t t_start1 = clock();
458  //Solve with ARPACK
459  {
461  problem(n_element);
462 
463  std::cout << "Matrix size " << problem.ndof() << std::endl;
464 
465  problem.solve(1);
466  }
467  clock_t t_end1 = clock();
468 
469  clock_t t_start2 = clock();
470  //Solve with LAPACK_QZ
471  {
473  problem(n_element);
474 
475  problem.solve(2);
476  }
477  clock_t t_end2 = clock();
478 
479 #ifdef OOMPH_HAS_TRILINOS
480  clock_t t_start3 = clock();
481 //Solve with Anasazi
482  {
483  HarmonicProblem<QHarmonicElement<3>,ANASAZI> problem(n_element);
484  problem.solve(3);
485  }
486  clock_t t_end3 = clock();
487 #endif
488 
489  std::cout << "ARPACK TIME: " << (double)(t_end1 - t_start1)/CLOCKS_PER_SEC
490  << std::endl;
491 
492  std::cout << "LAPACK TIME: " << (double)(t_end2 - t_start2)/CLOCKS_PER_SEC
493  << std::endl;
494 
495 #ifdef OOMPH_HAS_TRILINOS
496  std::cout << "ANASAZI TIME: " << (double)(t_end3 - t_start3)/CLOCKS_PER_SEC
497  << std::endl;
498 #endif
499 
500 #ifdef OOMPH_HAS_MPI
501  MPI_Helpers::finalize();
502 #endif
503 
504 } // end of main
505 
506 
507 
508 
509 
510 
511 
512 
513 
514 
unsigned required_nvalue(const unsigned &n) const
Required # of `values' (pinned or dofs) at node n.
Definition: harmonic.cc:238
double interpolated_u(const Vector< double > &s) const
Return FE representation of function value u(s) at local coordinate s.
Definition: harmonic.cc:174
QHarmonicElement()
Constructor: Call constructors for QElement and Poisson equations.
Definition: harmonic.cc:234
void output(ostream &outfile, const unsigned &Nplot)
Output function overloaded from HarmonicEquations.
Definition: harmonic.cc:245
HarmonicProblem(const unsigned &n_element)
Constructor: Pass number of elements and pointer to source function.
Definition: harmonic.cc:299
virtual double u(const unsigned &n) const
Access function: Eigenfunction value at local node n Note that solving the eigenproblem does not assi...
Definition: harmonic.cc:84
void output(ostream &outfile)
Output function overloaded from HarmonicEquations.
Definition: harmonic.cc:241
void solve(const unsigned &label)
Solve the problem.
Definition: harmonic.cc:361
int main(int argc, char **argv)
Driver for 1D Poisson problem.
Definition: harmonic.cc:446
bool operator()(const complex< T > &x, const complex< T > &y) const
Comparison. Are the values identical or not?
Definition: harmonic.cc:55
double dshape_eulerian_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsidx) const
Shape, test functions & derivs. w.r.t. to global coords. at integration point ipt. Return Jacobian.
Definition: harmonic.cc:260
void doc_solution(const unsigned &label)
Doc the solution, pass the number of the case considered, so that output files can be distinguished...
Definition: harmonic.cc:338
double dshape_eulerian(const Vector< double > &s, Shape &psi, DShape &dpsidx) const
Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
Definition: harmonic.cc:252
void output(ostream &outfile, const unsigned &nplot)
Output FE representation of soln: x,y,u or x,y,z,u at Nplot plot points.
Definition: harmonic.cc:96
virtual int u_local_eqn(const unsigned &n)
Access function that returns the local equation number of the unknown in the problem. Default is to assume that it is the first (only) value stored at the node.
Definition: harmonic.cc:211
void output(ostream &outfile)
Output the eigenfunction with default number of plot points.
Definition: harmonic.cc:88
1D Harmonic problem in unit interval.
Definition: harmonic.cc:272
void fill_in_contribution_to_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
Assemble the contributions to the jacobian and mass matrices.
Definition: harmonic.cc:119
~HarmonicProblem()
Destructor (empty)
Definition: harmonic.cc:280
HarmonicEquations()
Empty Constructor.
Definition: harmonic.cc:79