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