CppNoddy  0.92
Loading...
Searching...
No Matches
Classes | Namespaces | Macros | Functions
EVP2DOrrSommerfeld_slepcz_mumps.cpp File Reference
#include <cassert>
#include <Utility.h>
#include <SparseMatrix.h>
#include <SparseLinearEigenSystem.h>
#include <TwoD_Mapped_Node_Mesh.h>
#include <TwoD_Node_Mesh.h>
#include <SlepcSession.h>
#include <TrackerFile.h>

Go to the source code of this file.

Classes

class  Example::My_Mesh< _Type >
 

Namespaces

namespace  Example
 

Macros

#define UNIFORM
 
#define TYPEI
 
#define V_EVEN_Y
 
#define V_EVEN_Z
 
#define W_ODD_Y
 
#define W_ODD_Z
 

Functions

int main (int argc, char *argv[])
 

Macro Definition Documentation

◆ TYPEI

#define TYPEI

Definition at line 24 of file EVP2DOrrSommerfeld_slepcz_mumps.cpp.

◆ UNIFORM

#define UNIFORM

Definition at line 23 of file EVP2DOrrSommerfeld_slepcz_mumps.cpp.

◆ V_EVEN_Y

#define V_EVEN_Y

Definition at line 27 of file EVP2DOrrSommerfeld_slepcz_mumps.cpp.

◆ V_EVEN_Z

#define V_EVEN_Z

Definition at line 28 of file EVP2DOrrSommerfeld_slepcz_mumps.cpp.

◆ W_ODD_Y

#define W_ODD_Y

Definition at line 46 of file EVP2DOrrSommerfeld_slepcz_mumps.cpp.

◆ W_ODD_Z

#define W_ODD_Z

Definition at line 50 of file EVP2DOrrSommerfeld_slepcz_mumps.cpp.

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)

Definition at line 139 of file EVP2DOrrSommerfeld_slepcz_mumps.cpp.

140{
141 SlepcSession::getInstance(argc,argv);
142
143 cout << "\n";
144 cout << "=== EVP: Temporal spectra of the 2D Orr-Sommerfeld eqn =\n";
145 cout << "=== with a matrix problem assembled by hand using =\n";
146 cout << "=== the formulation of Tatsumi and Yoshimura (1990). =\n";
147 cout << "\n";
148
149 enum {v,w,F,G};
150
151 // duct aspect ratio
152 const double A(6.0);
153 // streamwise wavenumber and Reynolds number
154 const double alpha ( 0.94 );
155 double Re ( 8200.0 );
156 const D_complex eye( 0.0, 1.0 );
157
158 // discretise with these many nodal points
159 const std::size_t NZ( 301 );
160 const std::size_t NY( 301 );
161 //
162 //
163 Example::My_Mesh<double> base( -A, 0.0, -1.0, 0.0, NZ, NY, 1 );
164 //
165 for ( std::size_t i = 0; i < NZ; ++i )
166 {
167 for ( std::size_t j = 0; j < NY; ++j )
168 {
169 double Z = base.coord(i,j).first;
170 double Y = base.coord(i,j).second;
171 double sum = 1.0 - Y*Y;
172 for ( std::size_t n = 0; n < 1000; ++n )
173 {
174 double C = (2*n+1)*M_PI/2.0;
175 double term = cos(C*Y)*(exp(C*(Z-A))+exp(-C*(Z+A)))/(1+exp(-2*C*A));
176 sum -= 4*std::pow(2./M_PI,3)*std::pow(-1.0,n)*term/std::pow(2*n+1,3);
177 }
178 base(i,j,0) = sum;
179 }
180 }
181 base.dump_gnu("./DATA/U.dat");
182
183
184 // we'll solve as FOUR equations
185 const std::size_t N( 4*NZ*NY );
186
187 cout << "There are " << N << " degrees of freedom in this problem.\n";
188
189 // spatial step for the uniform computational-domain mesh
190 const double DZ = base.get_comp_step_sizes().first;
191 const double DY = base.get_comp_step_sizes().second;
192 //std::cout << " DZ = " << DZ << " DY = " << DY << "\n";
193
194 // matrices for the EVP, initialised with zeroes
195#ifdef DENSE
196 DenseMatrix<D_complex> a( N, N, 0.0 );
197 DenseMatrix<D_complex> b( N, N, 0.0 );
198#else
199 SparseMatrix<D_complex> a( N, N );
200 SparseMatrix<D_complex> b( N, N );
201#endif
202
203 std::size_t dof_per_i( 4*NY );
204 //
205 for ( std::size_t j = 0; j < NY; ++j )
206 {
207 const std::size_t i(0);
208 // left boundary conditions for the duct
209 const std::size_t O( i*4*NY );
210 std::size_t k( j*4 );
211 // coordinate mapping
212 double Z( base.coord(i,j).first );
213 double KZd( base.FnComp_Xd( Z ) );
214 double KZdd( base.FnComp_Xdd( Z ) );
215 // v = 0
216 a( O+k+0, O+k + v ) = 1.0;
217 // w = 0
218 a( O+k+1, O+k + w ) = 1.0;
219 // definition of F with w=0 for all y and 2nd order 1-sided diff.
220 a( O+k+2, O+k + dof_per_i + v ) = -5.0*(eye/(alpha*Re))*KZd*KZd/(DZ*DZ)
221 + 4.0*(eye/(alpha*Re))*KZdd/(2.0*DZ);
222 a( O+k+2, O+k + 2*dof_per_i + v ) = 4.0*(eye/(alpha*Re))*KZd*KZd/(DZ*DZ)
223 - 1.0*(eye/(alpha*Re))*KZdd/(2.0*DZ);
224 a( O+k+2, O+k + 3*dof_per_i + v ) = -1.0*(eye/(alpha*Re))*KZd*KZd/(DZ*DZ);
225 a( O+k+2, O+k + F ) = -1.0;
226 // definition of G -- building in w_z = 0 with high order scheme
227 a( O+k+3, O+k + dof_per_i + w ) = 10.0*(eye/(alpha*Re))*KZd*KZd/(3*DZ*DZ);
228 a( O+k+3, O+k + 2*dof_per_i + w ) = -1.0*(eye/(alpha*Re))*KZd*KZd/(3*DZ*DZ);
229 a( O+k+3, O+k + G ) = -1.0;
230 }
231 //
232 for ( std::size_t i = 1; i <= NZ - 2; ++i )
233 {
234 {
235 const std::size_t j(0);
236 // bottom boundary conditions for the duct
237 std::size_t O( i*4*NY );
238 std::size_t k( j*4 );
239 // coordinate mapping
240 double Y( base.coord(i,j).second );
241 double KYd( base.FnComp_Yd( Y ) );
242 double KYdd( base.FnComp_Ydd( Y ) );
243 // v=0 on y=0
244 a( O+k+0, O+k + v ) = 1.0;
245 // w=0 on y=0
246 a( O+k+1, O+k + w ) = 1.0;
247 // definition of F, building in v_y=0 with high order scheme
248 a( O+k+2, O+k + 4 + v ) = 10.0*(eye/(alpha*Re))*KYd*KYd/(3*DY*DY);
249 a( O+k+2, O+k + 8 + v ) = -1.0*(eye/(alpha*Re))*KYd*KYd/(3*DY*DY);
250 a( O+k+2, O+k + F ) = -1.0;
251 // definition of G with w=0 for all y and 2nd order 1-sided diff.
252 a( O+k+3, O+k + 4 + w ) = -5.0*(eye/(alpha*Re))*KYd*KYd/(DY*DY)
253 + 4.0*(eye/(alpha*Re))*KYdd/(2.0*DY);
254 a( O+k+3, O+k + 8 + w ) = 4.0*(eye/(alpha*Re))*KYd*KYd/(DY*DY)
255 - 1.0*(eye/(alpha*Re))*KYdd/(2.0*DY);
256 a( O+k+3, O+k + 12 + w ) = -1.0*(eye/(alpha*Re))*KYd*KYd/(DY*DY);
257 a( O+k+3, O+k + G ) = -1.0;
258 }
259 for ( std::size_t j = 1; j < NY-1; ++j )
260 {
261 // coordinate mapping
262 double Z( base.coord(i,j).first );
263 double KZd( base.FnComp_Xd( Z ) );
264 double KZdd( base.FnComp_Xdd( Z ) );
265 double Y( base.coord(i,j).second );
266 double KYd( base.FnComp_Yd( Y ) );
267 double KYdd( base.FnComp_Ydd( Y ) );
268 //
269 double Uy = (base(i,j+1,0)-base(i,j-1,0))*KYd/(2*DY);
270 double Uz = (base(i+1,j,0)-base(i-1,j,0))*KZd/(2*DZ);
271 double Uyy = (base(i,j+1,0)-2*base(i,j,0)+base(i,j-1,0))*KYd*KYd/(DY*DY)
272 + (base(i,j+1,0)-base(i,j-1,0))*KYdd/(2*DY);
273 double Uzz = (base(i+1,j,0)-2*base(i,j,0)+base(i-1,j,0))*KZd*KZd/(DZ*DZ)
274 + (base(i+1,j,0)-base(i-1,j,0))*KZdd/(2*DZ);
275 double Uyz = (base(i+1,j+1,0)-base(i+1,j-1,0)
276 - base(i-1,j+1,0)+base(i-1,j-1,0))*KYd*KZd/(4*DY*DZ);
277 // interior equations within the duct
278 std::size_t O( i*4*NY );
279 std::size_t k( j*4 );
280 //
281 // L{v}=F
282 //
283 // v_i,j
284 a( O+k+0, O+k + v ) = (eye/(alpha*Re))*( - 2.0*KYd*KYd/(DY*DY)
285 - 2.0*KZd*KZd/(DZ*DZ)
286 - alpha*alpha ) + base(i,j,0);
287 // v_i+1,j
288 a( O+k+0, O+k + dof_per_i + v ) = (eye/(alpha*Re))*KZd*KZd/(DZ*DZ)
289 + (eye/(alpha*Re))*KZdd/(2.0*DZ);
290 // v_i-1,j
291 a( O+k+0, O+k - dof_per_i + v ) = (eye/(alpha*Re))*KZd*KZd/(DZ*DZ)
292 - (eye/(alpha*Re))*KZdd/(2.0*DZ);
293 // v_i,j+1
294 a( O+k+0, O+k + 4 + v ) = (eye/(alpha*Re))*KYd*KYd/(DY*DY)
295 + (eye/(alpha*Re))*KYdd/(2.0*DY);
296 // v_i,j-1
297 a( O+k+0, O+k - 4 + v ) = (eye/(alpha*Re))*KYd*KYd/(DY*DY)
298 - (eye/(alpha*Re))*KYdd/(2.0*DY);
299 // F_i,j
300 a( O+k+0, O+k + F ) = -1.0;
301 // ev -- wave speed c
302 b( O+k+0, O+k + v ) = 1.0;
303 //
304 // L{w}=G
305 //
306 // w_i,j
307 a( O+k+1, O+k + w ) = (eye/(alpha*Re))*( - 2.0*KYd*KYd/(DY*DY)
308 - 2.0*KZd*KZd/(DZ*DZ)
309 - alpha*alpha ) + base(i,j,0);
310 // w_i+1,j
311 a( O+k+1, O+k + dof_per_i + w ) = (eye/(alpha*Re))*KZd*KZd/(DZ*DZ)
312 + (eye/(alpha*Re))*KZdd/(2.0*DZ);
313 // w_i-1,j
314 a( O+k+1, O+k - dof_per_i + w ) = (eye/(alpha*Re))*KZd*KZd/(DZ*DZ)
315 - (eye/(alpha*Re))*KZdd/(2.0*DZ);
316 // w_i,j+1
317 a( O+k+1, O+k + 4 + w ) = (eye/(alpha*Re))*KYd*KYd/(DY*DY)
318 + (eye/(alpha*Re))*KYdd/(2.0*DY);
319 // w_i,j-1
320 a( O+k+1, O+k - 4 + w ) = (eye/(alpha*Re))*KYd*KYd/(DY*DY)
321 - (eye/(alpha*Re))*KYdd/(2.0*DY);
322 // G_i,j
323 a( O+k+1, O+k + G ) = -1.0;
324 // ev -- wave speed c
325 b( O+k+1, O+k + w ) = 1.0;
326 //
327 // E(y,z){v} = O(y,z){w}
328 //
329 // F_i,j
330 a( O+k+2, O+k + F ) = -2.0*KYd*KYd/(DY*DY) - alpha*alpha;
331 // F_i,j+1
332 a( O+k+2, O+k + 4 + F ) = 1.0*KYd*KYd/(DY*DY) + KYdd/(2.0*DY);
333 // F_i,j-1
334 a( O+k+2, O+k - 4 + F ) = 1.0*KYd*KYd/(DY*DY) - KYdd/(2.0*DY);
335 //
336 // G_i+1,j+1
337 a( O+k+2, O+k + dof_per_i + 4 + G ) = 0.25*KYd*KZd/(DY*DZ);
338 // G_i+1,j-1
339 a( O+k+2, O+k + dof_per_i - 4 + G ) = -0.25*KYd*KZd/(DY*DZ);
340 // G_i-1,j-1
341 a( O+k+2, O+k - dof_per_i - 4 + G ) = 0.25*KYd*KZd/(DY*DZ);
342 // G_i-1,j+1
343 a( O+k+2, O+k - dof_per_i + 4 + G ) = -0.25*KYd*KZd/(DY*DZ);
344 //
345 // v_i,j
346 a( O+k+2, O+k + v ) = -2.0*Uyy;
347 // v_i,j+1
348 a( O+k+2, O+k + 4 + v ) = -2.0*Uy*KYd/(2.0*DY);
349 // v_i,j-1
350 a( O+k+2, O+k - 4 + v ) = +2.0*Uy*KYd/(2.0*DY);
351 //
352 // w_i,j
353 a( O+k+2, O+k + w ) = -2.0*Uyz;
354 // w_i,j+1
355 a( O+k+2, O+k + 4 + w ) = -2.0*Uz*KYd/(2.0*DY);
356 // w_i,j-1
357 a( O+k+2, O+k - 4 + w ) = +2.0*Uz*KYd/(2.0*DY);
358 //
359 //
360 // E(z,y){w} = O(z,y){v}
361 //
362 // G_i,j
363 a( O+k+3, O+k + G ) = -2.0*KZd*KZd/(DZ*DZ) - alpha*alpha;
364 // G_i+1,j
365 a( O+k+3, O+k + dof_per_i + G ) = 1.0*KZd*KZd/(DZ*DZ) + KZdd/(2.0*DZ);
366 // G_i-1,j
367 a( O+k+3, O+k - dof_per_i + G ) = 1.0*KZd*KZd/(DZ*DZ) - KZdd/(2.0*DZ);
368 //
369 // F_i+1,j+1
370 a( O+k+3, O+k + dof_per_i + 4 + F ) = 0.25*KYd*KZd/(DY*DZ);
371 // F_i+1,j-1
372 a( O+k+3, O+k + dof_per_i - 4 + F ) = -0.25*KYd*KZd/(DY*DZ);
373 // F_i-1,j-1
374 a( O+k+3, O+k - dof_per_i - 4 + F ) = 0.25*KYd*KZd/(DY*DZ);
375 // F_i-1,j+1
376 a( O+k+3, O+k - dof_per_i + 4 + F ) = -0.25*KYd*KZd/(DY*DZ);
377 //
378 // v_i,j
379 a( O+k+3, O+k + v ) = -2.0*Uyz;
380 // v_i+1,j
381 a( O+k+3, O+k + dof_per_i + v ) = -2.0*Uy*KZd/(2.0*DZ);
382 // v_i-1,j
383 a( O+k+3, O+k - dof_per_i + v ) = +2.0*Uy*KZd/(2.0*DZ);
384 //
385 // w_i,j
386 a( O+k+3, O+k + w ) = -2.0*Uzz;
387 // w_i+1,j
388 a( O+k+3, O+k + dof_per_i + w ) = -2.0*Uz*KZd/(2.0*DZ);
389 // w_i-1,j
390 a( O+k+3, O+k - dof_per_i + w ) = +2.0*Uz*KZd/(2.0*DZ);
391 }
392 {
393 // top boundary
394 //
395 const std::size_t j(NY-1);
396 std::size_t O( i*4*NY );
397 const std::size_t k( j*4 );
398 // coordinate mapping
399 double Y( base.coord(i,j).second );
400 double KYd( base.FnComp_Yd( Y ) );
401 //double KYdd( base.FnComp_Ydd( Y ) );
402 //
403#ifdef V_EVEN_Y
404 // v_y = 0
405 a( O+k+0, O+k + v ) = 3.0*KYd/(2.0*DY);
406 a( O+k+0, O+k - 4 + v ) = -4.0*KYd/(2.0*DY);
407 a( O+k+0, O+k - 8 + v ) = 1.0*KYd/(2.0*DY);
408 // v_yyy = 0 => F_y = 0
409 a( O+k+2, O+k + F ) = 3.0*KYd/(2.0*DY);
410 a( O+k+2, O+k - 4 + F ) = -4.0*KYd/(2.0*DY);
411 a( O+k+2, O+k - 8 + F ) = 1.0*KYd/(2.0*DY);
412#elif defined(V_ODD_Y)
413 // v=0
414 a( O+k+0, O+k + v ) = 1.0;
415 // v_yy =0 => F=0
416 a( O+k+2, O+k + F ) = 1.0;
417#endif
418
419#ifdef W_EVEN_Y
420 // w_y = 0
421 a( O+k+1, O+k + w ) = 3.0*KYd/(2.0*DY);
422 a( O+k+1, O+k - 4 + w ) = -4.0*KYd/(2.0*DY);
423 a( O+k+1, O+k - 8 + w ) = 1.0*KYd/(2.0*DY);
424 // w_yyy = 0 => G_y = 0
425 a( O+k+3, O+k + G ) = 3.0*KYd/(2.0*DY);
426 a( O+k+3, O+k - 4 + G ) = -4.0*KYd/(2.0*DY);
427 a( O+k+3, O+k - 8 + G ) = 1.0*KYd/(2.0*DY);
428#elif defined(W_ODD_Y)
429 // w = 0
430 a( O+k+1, O+k + w ) = 1.0;
431 // G = 0
432 a( O+k+3, O+k + G ) = 1.0;
433#endif
434 }
435 }
436 {
437 for ( std::size_t j = 0; j < NY; ++j )
438 {
439 // right hand boundary
440 const std::size_t i(NZ-1);
441 // left boundary conditions for the duct
442 const std::size_t O( i*4*NY );
443 std::size_t k( j*4 );
444 // coordinate mapping
445 double Z( base.coord(i,j).first );
446 double KZd( base.FnComp_Xd( Z ) );
447 //double KZdd( base.FnComp_Xdd( Z ) );
448#ifdef V_EVEN_Z
449 // v_z = 0
450 a( O+k+0, O+k + v ) = 3.0*KZd/(2.0*DZ);
451 a( O+k+0, O+k - dof_per_i + v ) = -4.0*KZd/(2.0*DZ);
452 a( O+k+0, O+k - 2*dof_per_i + v ) = 1.0*KZd/(2.0*DZ);
453 // v_zzz = 0 => F_z = 0
454 a( O+k+2, O+k + F ) = 3.0*KZd/(2.0*DZ);
455 a( O+k+2, O+k - dof_per_i + F ) = -4.0*KZd/(2.0*DZ);
456 a( O+k+2, O+k - 2*dof_per_i + F ) = 1.0*KZd/(2.0*DZ);
457#elif defined(V_ODD_Z)
458 // v=0
459 a( O+k+0, O+k + v ) = 1.0;
460 // v_zz =0 => F=0
461 a( O+k+2, O+k + F ) = 1.0;
462#endif
463
464#ifdef W_EVEN_Z
465 // w_z = 0
466 a( O+k+1, O+k + w ) = 3.0*KZd/(2.0*DZ);
467 a( O+k+1, O+k - dof_per_i + w ) = -4.0*KZd/(2.0*DZ);
468 a( O+k+1, O+k - 2*dof_per_i + w ) = 1.0*KZd/(2.0*DZ);
469 // w_zzz = 0 => G_z = 0
470 a( O+k+3, O+k + G ) = 3.0*KZd/(2.0*DZ);
471 a( O+k+3, O+k - dof_per_i + G ) = -4.0*KZd/(2.0*DZ);
472 a( O+k+3, O+k - 2*dof_per_i + G ) = 1.0*KZd/(2.0*DZ);
473#elif defined(W_ODD_Z)
474 // w = 0
475 a( O+k+1, O+k + w ) = 1.0;
476 // w_zz = 0 => G = 0
477 a( O+k+3, O+k + G ) = 1.0;
478#endif
479 }
480 }
481
482 // a vector for the eigenvalues
484#ifdef DENSE
485 DenseLinearEigenSystem<D_complex> system( &a, &b );
486#else
487 SparseLinearEigenSystem<D_complex> system( &a, &b );
488 system.set_nev(1);
489 system.set_region(0.0, 1.0, -1.0, 1.0);
490 system.set_target( D_complex( 0.244, 0.0 ) );
491 system.set_order("EPS_TARGET_IMAGINARY");
492#endif
493
494 system.eigensolve();
495
496 system.tag_eigenvalues_disc( + 1, 1.0 );
497 ev_c = system.get_tagged_eigenvalues();
498 TrackerFile spectrum( "./DATA/spectrum_disc_2d.dat" );
499 spectrum.push_ptr( &ev_c, "evs" );
500 spectrum.update();
501 ev_c.dump();
502
503 DenseMatrix<D_complex> vec_mtx = system.get_tagged_eigenvectors();
504 Example::My_Mesh<D_complex> eigfn( -A, 0.0, -1.0, 0.0, NZ, NY, 4 );
505
506 for ( std::size_t n = 0; n < ev_c.size(); ++n )
507 {
508 for ( std::size_t i = 0; i < NZ; ++i )
509 {
510 for ( std::size_t j = 0; j < NY; ++j )
511 {
512 eigfn(i,j,v) = vec_mtx(n, i*4*NY + 4*j + v );
513 eigfn(i,j,w) = vec_mtx(n, i*4*NY + 4*j + w );
514 eigfn(i,j,F) = vec_mtx(n, i*4*NY + 4*j + F );
515 eigfn(i,j,G) = vec_mtx(n, i*4*NY + 4*j + G );
516 }
517 }
518 eigfn.dump_gnu("./DATA/eigfn_"+Utility::stringify(n,4)+".dat");
519 }
520
521 if ( std::abs( ev_c[0].imag() ) < 1.e-4 )
522 {
523 cout << "\033[1;32;48m * PASSED \033[0m\n";
524 return 0;
525 }
526
527 cout << "\033[1;31;48m * FAILED \033[0m\n";
528 cout << " Final error = " << ev_c[0].imag() << "\n";
529 return 1;
530
531}
A linear Nth-order generalised eigensystem class.
A matrix class that constructs a DENSE matrix as a row major std::vector of DenseVectors.
Definition: DenseMatrix.h:25
An DenseVector class – a dense vector object.
Definition: DenseVector.h:34
void dump() const
Dump to std::cout.
Definition: DenseVector.cpp:64
std::size_t size() const
A pass-thru definition to get the size of the vector.
Definition: DenseVector.h:330
A matrix class that constructs a SPARSE matrix as a row major std::vector of SparseVectors.
Definition: SparseMatrix.h:31
double A(1.0)
initial hump amplitude
const D_complex eye(0., 1.0)
double Re
Globally define the Reynolds number and wavenumber.
DenseVector< double > imag(const DenseVector< D_complex > &X)
Return a double DENSE vector containing the imaginary part of a complex DENSE vector.
Definition: Utility.cpp:185
std::string stringify(const int &val)
Return an integer value as a string - useful for file naming.
Definition: Utility.cpp:193
std::complex< double > D_complex
A complex double precision number using std::complex.
Definition: Types.h:98

References CppNoddy::TwoD_Mapped_Node_Mesh< _Type >::coord(), CppNoddy::DenseVector< _Type >::dump(), CppNoddy::TwoD_Mapped_Node_Mesh< _Type >::dump_gnu(), CppNoddy::DenseLinearEigenSystem< _Type >::eigensolve(), Example::My_Mesh< _Type >::FnComp_Xd(), Example::My_Mesh< _Type >::FnComp_Xdd(), Example::My_Mesh< _Type >::FnComp_Yd(), Example::My_Mesh< _Type >::FnComp_Ydd(), CppNoddy::TwoD_Mapped_Node_Mesh< _Type >::get_comp_step_sizes(), CppNoddy::DenseLinearEigenSystem< _Type >::get_tagged_eigenvalues(), CppNoddy::DenseLinearEigenSystem< _Type >::get_tagged_eigenvectors(), CppNoddy::TrackerFile::push_ptr(), CppNoddy::DenseVector< _Type >::size(), CppNoddy::Utility::stringify(), CppNoddy::DenseLinearEigenSystem< _Type >::tag_eigenvalues_disc(), and CppNoddy::TrackerFile::update().

© 2012

R.E. Hewitt