From: David Xu <dxu@my...>  20060304 18:14:50

Hi Steffen, I'm not familiar with Slepc. Could elaborate how to modify the code or use command line options to calculate the smallest eigenvalues? Thanks, David On 2/27/06, Steffen Petersen <steffen.petersen@...> wrote: > I have just committed an example 17 including a generalized > eigenvalue problem. Currently, the largest eigenvalues > in the spectrum are computed. The code can be easily changed > in order to compute the smallest values (or use the slepc comman > line options). However, the Lanczos solver slepc offers with this > option appears to be somewhat unstable when used to compute > the smallest values (perhaps this is due to the zero eigenmode in > this example). An option to work around this problem is to > configure and install slepc with arpack and parpack (the arpack > solver works fine for example 17). > > Steffen > > > > David Xu schrieb: > > Hi Steffen and Ondrej, > > > > I read your previous message that you've extended the slepc interface > > to support generalized eigen problems. I'd like to get the details > > from you. I'm new to libmesh and I'm wondering if you can post some > > new tutorial examples on how to solve eigen problems using libmesh > > (eg. solve 1D or 2D simple harmonic oscillator)? That will a great > > help. > > > > Ondrej, I read about your message about your implementation of libmesh > > and pysparse's jacobidavidson solver. That sounds really good. Could > > you please give me the details too? I really appreciate it. > > > > David > > 
From: David Xu <dxu@my...>  20060220 22:48:43

Hi Steffen and Ondrej, I read your previous message that you've extended the slepc interface to support generalized eigen problems. I'd like to get the details from you. I'm new to libmesh and I'm wondering if you can post some new tutorial examples on how to solve eigen problems using libmesh (eg. solve 1D or 2D simple harmonic oscillator)? That will a great help. Ondrej, I read about your message about your implementation of libmesh and pysparse's jacobidavidson solver. That sounds really good. Could you please give me the details too? I really appreciate it. David 
From: Steffen Petersen <steffen.petersen@tu...>  20060221 18:28:34

Hello David, I will commit an adjusted example within the next days. Note that some rather useful functionality is only available with the latest slepc version (e.g. computation of the smallest eigenvalues within the spectrum), so I recommend to use the 2.3.0 versions of petsc and slepc. Steffen > Hi Steffen and Ondrej, > > I read your previous message that you've extended the slepc interface > to support generalized eigen problems. I'd like to get the details > from you. I'm new to libmesh and I'm wondering if you can post some > new tutorial examples on how to solve eigen problems using libmesh > (eg. solve 1D or 2D simple harmonic oscillator)? That will a great > help. > > Ondrej, I read about your message about your implementation of libmesh > and pysparse's jacobidavidson solver. That sounds really good. Could > you please give me the details too? I really appreciate it. > > David > 
From: Steffen Petersen <steffen.petersen@tu...>  20060227 08:33:35

I have just committed an example 17 including a generalized eigenvalue problem. Currently, the largest eigenvalues in the spectrum are computed. The code can be easily changed in order to compute the smallest values (or use the slepc comman line options). However, the Lanczos solver slepc offers with this option appears to be somewhat unstable when used to compute the smallest values (perhaps this is due to the zero eigenmode in this example). An option to work around this problem is to configure and install slepc with arpack and parpack (the arpack solver works fine for example 17). Steffen David Xu schrieb: > Hi Steffen and Ondrej, > > I read your previous message that you've extended the slepc interface > to support generalized eigen problems. I'd like to get the details > from you. I'm new to libmesh and I'm wondering if you can post some > new tutorial examples on how to solve eigen problems using libmesh > (eg. solve 1D or 2D simple harmonic oscillator)? That will a great > help. > > Ondrej, I read about your message about your implementation of libmesh > and pysparse's jacobidavidson solver. That sounds really good. Could > you please give me the details too? I really appreciate it. > > David 
From: David Xu <dxu@my...>  20060301 18:59:50

Hi Steffen, I'm looking forward to your example 17 and implementation of generalized eigen value solver with libmesh. Thanks, David On 2/27/06, Steffen Petersen <steffen.petersen@...> wrote: > I have just committed an example 17 including a generalized > eigenvalue problem. Currently, the largest eigenvalues > in the spectrum are computed. The code can be easily changed > in order to compute the smallest values (or use the slepc comman > line options). However, the Lanczos solver slepc offers with this > option appears to be somewhat unstable when used to compute > the smallest values (perhaps this is due to the zero eigenmode in > this example). An option to work around this problem is to > configure and install slepc with arpack and parpack (the arpack > solver works fine for example 17). > > Steffen > > > > David Xu schrieb: > > Hi Steffen and Ondrej, > > > > I read your previous message that you've extended the slepc interface > > to support generalized eigen problems. I'd like to get the details > > from you. I'm new to libmesh and I'm wondering if you can post some > > new tutorial examples on how to solve eigen problems using libmesh > > (eg. solve 1D or 2D simple harmonic oscillator)? That will a great > > help. > > > > Ondrej, I read about your message about your implementation of libmesh > > and pysparse's jacobidavidson solver. That sounds really good. Could > > you please give me the details too? I really appreciate it. > > > > David > > 
From: David Xu <dxu@my...>  20060304 18:14:50

Hi Steffen, I'm not familiar with Slepc. Could elaborate how to modify the code or use command line options to calculate the smallest eigenvalues? Thanks, David On 2/27/06, Steffen Petersen <steffen.petersen@...> wrote: > I have just committed an example 17 including a generalized > eigenvalue problem. Currently, the largest eigenvalues > in the spectrum are computed. The code can be easily changed > in order to compute the smallest values (or use the slepc comman > line options). However, the Lanczos solver slepc offers with this > option appears to be somewhat unstable when used to compute > the smallest values (perhaps this is due to the zero eigenmode in > this example). An option to work around this problem is to > configure and install slepc with arpack and parpack (the arpack > solver works fine for example 17). > > Steffen > > > > David Xu schrieb: > > Hi Steffen and Ondrej, > > > > I read your previous message that you've extended the slepc interface > > to support generalized eigen problems. I'd like to get the details > > from you. I'm new to libmesh and I'm wondering if you can post some > > new tutorial examples on how to solve eigen problems using libmesh > > (eg. solve 1D or 2D simple harmonic oscillator)? That will a great > > help. > > > > Ondrej, I read about your message about your implementation of libmesh > > and pysparse's jacobidavidson solver. That sounds really good. Could > > you please give me the details too? I really appreciate it. > > > > David > > 
From: David Xu <dxu@my...>  20060306 08:17:22

Hi Steffen, Besides ARPACK, can I call other package (BLZPACK, PLANSO, and TRLAN) from libmesh that SLEPC has an interface to? Of course I will install them with SLEPC first. Thanks, David On 2/27/06, Steffen Petersen <steffen.petersen@...> wrote: > I have just committed an example 17 including a generalized > eigenvalue problem. Currently, the largest eigenvalues > in the spectrum are computed. The code can be easily changed > in order to compute the smallest values (or use the slepc comman > line options). However, the Lanczos solver slepc offers with this > option appears to be somewhat unstable when used to compute > the smallest values (perhaps this is due to the zero eigenmode in > this example). An option to work around this problem is to > configure and install slepc with arpack and parpack (the arpack > solver works fine for example 17). > > Steffen > > > > David Xu schrieb: > > Hi Steffen and Ondrej, > > > > I read your previous message that you've extended the slepc interface > > to support generalized eigen problems. I'd like to get the details > > from you. I'm new to libmesh and I'm wondering if you can post some > > new tutorial examples on how to solve eigen problems using libmesh > > (eg. solve 1D or 2D simple harmonic oscillator)? That will a great > > help. > > > > Ondrej, I read about your message about your implementation of libmesh > > and pysparse's jacobidavidson solver. That sounds really good. Could > > you please give me the details too? I really appreciate it. > > > > David > > 
From: Steffen Petersen <steffen.petersen@tu...>  20060302 10:26:06

Hello David, example 17 is already available. If you download the latest cvs version with cvs d:pserver:anonymous@...:/cvsroot/libmesh co libmesh you will find ex17 in the example directory. If you want to update your current version, you will have to update libmesh with the "d" option cvs d:pserver:ano... update d libmesh in order to receive the files that have been added recently. The online documentation does not include ex17 since it is not included in the actual libmesh release (although I think most people use the cvs version). As far as I know, the online documentation corresponds to the latest libmesh release (someone may correct me if I am wrong). Anyway, I will commit the corresponding files for the web page soon. Steffen David Xu schrieb: > Hi Steffen, > > I'm looking forward to your example 17 and implementation of > generalized eigen value solver with libmesh. > > Thanks, > > David > > On 2/27/06, Steffen Petersen <steffen.petersen@...> wrote: > >>I have just committed an example 17 including a generalized >>eigenvalue problem. Currently, the largest eigenvalues >>in the spectrum are computed. The code can be easily changed >>in order to compute the smallest values (or use the slepc comman >>line options). However, the Lanczos solver slepc offers with this >>option appears to be somewhat unstable when used to compute >>the smallest values (perhaps this is due to the zero eigenmode in >>this example). An option to work around this problem is to >>configure and install slepc with arpack and parpack (the arpack >>solver works fine for example 17). >> >>Steffen >> >> >> >>David Xu schrieb: >> >>>Hi Steffen and Ondrej, >>> >>>I read your previous message that you've extended the slepc interface >>>to support generalized eigen problems. I'd like to get the details >>>from you. I'm new to libmesh and I'm wondering if you can post some >>>new tutorial examples on how to solve eigen problems using libmesh >>>(eg. solve 1D or 2D simple harmonic oscillator)? That will a great >>>help. >>> >>>Ondrej, I read about your message about your implementation of libmesh >>>and pysparse's jacobidavidson solver. That sounds really good. Could >>>you please give me the details too? I really appreciate it. >>> >>>David >> >> 
From: David Xu <dxu@my...>  20060303 09:07:52

Steffen, Thanks for the heads up. I found the ex17 in the cvs. I question is how to export the desired eigenvalues to a file or standard output? You have the following lines in ex16 and ex17, which is to output the eigenvectors to a gmv file: =09eigen_system.get_eigenpair(nconv1); =09 =09// Write the eigen vector to file. =09char buf[14]; =09sprintf (buf, "out.gmv"); =09GMVIO (mesh).write_equation_systems (buf, equation_systems); I didn't how to use eigen_system.get_eigenpair to get the eigenvalues. Please enlighten me. Also, please let me know how to export the system matrices A and B to a file, Ondrej Certik used pysparse/JacobiDavidson program on the system matrices after they are exported to files. I'd like to try that too. Thank you, David On 3/2/06, Steffen Petersen <steffen.petersen@...> wrote: > > Hello David, > > example 17 is already available. If you download the latest > cvs version with > cvs d:pserver:anonymous@...:/cvsroot/libmesh co libmesh > you will find ex17 in the example directory. > If you want to update your current version, you will > have to update libmesh with the "d" option > cvs d:pserver:ano... update d libmesh > in order to receive the files that have been added recently. > > The online documentation does not include ex17 > since it is not included in the actual libmesh release > (although I think most people use the cvs version). > As far as I know, the online documentation corresponds to > the latest libmesh release (someone may correct me if I am wrong). > Anyway, I will commit the corresponding files for the > web page soon. > > Steffen > > > David Xu schrieb: > > Hi Steffen, > > > > I'm looking forward to your example 17 and implementation of > > generalized eigen value solver with libmesh. > > > > Thanks, > > > > David > > > > On 2/27/06, Steffen Petersen <steffen.petersen@...> wrote: > > > >>I have just committed an example 17 including a generalized > >>eigenvalue problem. Currently, the largest eigenvalues > >>in the spectrum are computed. The code can be easily changed > >>in order to compute the smallest values (or use the slepc comman > >>line options). However, the Lanczos solver slepc offers with this > >>option appears to be somewhat unstable when used to compute > >>the smallest values (perhaps this is due to the zero eigenmode in > >>this example). An option to work around this problem is to > >>configure and install slepc with arpack and parpack (the arpack > >>solver works fine for example 17). > >> > >>Steffen > >> > >> > >> > >>David Xu schrieb: > >> > >>>Hi Steffen and Ondrej, > >>> > >>>I read your previous message that you've extended the slepc interface > >>>to support generalized eigen problems. I'd like to get the details > >>>from you. I'm new to libmesh and I'm wondering if you can post some > >>>new tutorial examples on how to solve eigen problems using libmesh > >>>(eg. solve 1D or 2D simple harmonic oscillator)? That will a great > >>>help. > >>> > >>>Ondrej, I read about your message about your implementation of libmesh > >>>and pysparse's jacobidavidson solver. That sounds really good. Could > >>>you please give me the details too? I really appreciate it. > >>> > >>>David > >> > >> > > 
From: David Xu <dxu@my...>  20060304 09:59:38

Hi Everyone, I followed Steffen's new generalized eigen problem code in the new ex17 and wrote my first libmesh code to solve 1D timeindependent schrodinger equation with simple harmonic oscillator potential, the equation is 0.5u''+0.5*x^2*u=3DE*u . The analytical solution of the eigenvalue E are: 0.5, 1.5, 2.5, 3.5. 4.5 ... I set up the 1D domain [5 5] and apply the Dirichlet BC u(5)=3Du(5)=3D0 following ex0. The code compiled okay, but the results were way large eigenvalues and I don't know what went wrong in my code. So, I'm attaching my code here, it should fairly straightforward since it's a simple 1D case. I'd really appreciate it if anyone can help me figure out the problem. Source code: =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D // 1D Time independent Schrodinger Equation Eigenvalue Solver // using simple harmonic potential // libMesh include files. #include "libmesh.h" #include "mesh.h" #include "mesh_generation.h" #include "gmv_io.h" #include "eigen_system.h" #include "equation_systems.h" #include "fe.h" #include "quadrature_gauss.h" #include "dense_matrix.h" #include "sparse_matrix.h" #include "numeric_vector.h" #include "dof_map.h" #include "edge_edge3.h" // Function prototype. This is the function that will assemble // the eigen system. Here, we will simply assemble a mass matrix. void assemble_mass(EquationSystems& es, =09=09 const std::string& system_name); int main (int argc, char** argv) { // Initialize libMesh and the dependent libraries. libMesh::init (argc, argv); // This example is designed for the SLEPc eigen solver interface. #ifndef HAVE_SLEPC std::cerr << "ERROR: This example requires libMesh to be\n" =09 << "compiled with SLEPc eigen solvers support!" =09 << std::endl; return 0; #else // Braces are used to force object scope. { // Check for proper usage. if (argc < 3) { =09std::cerr << "\nUsage: " << argv[0] =09=09 << " n <number of eigen values>" =09=09 << std::endl; =09error(); } // Tell the user what we are doing. else { =09std::cout << "Running " << argv[0]; =09 =09for (int i=3D1; i<argc; i++) =09 std::cout << " " << argv[i]; =09 =09std::cout << std::endl << std::endl; } // Set the dimensionality. const unsigned int dim =3D 1; // Get the number of eigen values to be computed from argv[2] const unsigned int nev =3D std::atoi(argv[2]); // Create a dimdimensional mesh. Mesh mesh (dim); // Use the internal mesh generator to create a uniform // grid on a 1D [5 5] domain. MeshTools::Generation::build_line (mesh, 100, 5., 5., EDGE3); // Print information about the mesh to the screen. mesh.print_info(); // Create an equation systems object. EquationSystems equation_systems (mesh); // Create a EigenSystem named "Eigensystem" and (for convenience) // use a reference to the system we create. EigenSystem & eigen_system =3D equation_systems.add_system<EigenSystem> ("Eigensystem"); // Declare the system variables. { // Adds the variable "p" to "Eigensystem". "p" // will be approximated using secondorder approximation. eigen_system.add_variable("u", SECOND); // Give the system a pointer to the matrix assembly // function defined below. eigen_system.attach_assemble_function (assemble_mass); // Set necessary parametrs used in EigenSystem::solve(), // i.e. the number of requested eigenpairs \p nev and the number // of basis vectors \p ncv used in the solution algorithm. Note that // ncv >=3D nev must hold and ncv >=3D 2*nev is recommended. equation_systems.parameters.set<unsigned int>("eigenpairs") =3D ne= v; equation_systems.parameters.set<unsigned int>("basis vectors") =3D ne= v*3; // Set the eigen solver type. SLEPc offers various solvers such as // the Arnoldi and subspace method. It // also offers interfaces to other solver packages (e.g. ARPACK). eigen_system.eigen_solver>set_eigensolver_type(ARNOLDI); // Set the solver tolerance and the maximum number of iterations. equation_systems.parameters.set<Real>("linear solver tolerance") =3D pow(TOLERANCE, 5./3.); equation_systems.parameters.set<unsigned int> =09("linear solver maximum iterations") =3D 1000; // Set the type of the problem, here we deal with // a generalized Hermitian problem. eigen_system.set_eigenproblem_type(GHEP); // Set the eigenvalues to be computed. Note that not // all solvers support this. // eigen_system.eigen_solver>set_position_of_spectrum(SMALLEST_MAGNI= TUDE); // Initialize the data structures for the equation system. equation_systems.init(); // Prints information about the system to the screen. equation_systems.print_info(); } // Solve the system "Eigensystem". eigen_system.solve(); // Get the number of converged eigen pairs. unsigned int nconv =3D eigen_system.get_n_converged(); std::cout << "Number of converged eigenpairs: " << nconv =09 << "\n" << std::endl; std::pair<Real, Real> myeig; // Get the last converged eigenpair if (nconv !=3D 0) { =09//print eigenvalues for(unsigned int n=3D0; n<=3Dnconv1; n++) { =09myeig=3Deigen_system.get_eigenpair(n); std::cout << "eig(" << n << ")=3D " << myeig.first << "\n" << std::= endl; } } else { =09std::cout << "WARNING: Solver did not converge!\n" << nconv << std::endl= ; } } #endif // HAVE_SLEPC // All done. return libMesh::close (); } void assemble_mass(EquationSystems& es, =09=09 const std::string& system_name) { // It is a good idea to make sure we are assembling // the proper system. assert (system_name =3D=3D "Eigensystem"); #ifdef HAVE_SLEPC // Get a constant reference to the mesh object. const Mesh& mesh =3D es.get_mesh(); // The dimension that we are running. const unsigned int dim =3D mesh.mesh_dimension(); // Get a reference to our system. EigenSystem & eigen_system =3D es.get_system<EigenSystem> (system_name); // Get a constant reference to the Finite Element type // for the first (and only) variable in the system. FEType fe_type =3D eigen_system.get_dof_map().variable_type(0); // A reference to the two system matrices SparseMatrix<Number>& matrix_A =3D *eigen_system.matrix_A; SparseMatrix<Number>& matrix_B =3D *eigen_system.matrix_B; // Build a Finite Element object of the specified type. Since the // \p FEBase::build() member dynamically creates memory we will // store the object as an \p AutoPtr<FEBase>. This can be thought // of as a pointer that will clean up after itself. AutoPtr<FEBase> fe (FEBase::build(dim, fe_type)); // A Gauss quadrature rule for numerical integration. // Use the default quadrature order. //QGauss qrule (dim, fe_type.default_quadrature_order()); QGauss qrule(dim,FIFTH); // Tell the finite element object to use our quadrature rule. fe>attach_quadrature_rule (&qrule); // The element Jacobian * quadrature weight at each integration point. const std::vector<Real>& JxW =3D fe>get_JxW(); const std::vector<Point>& q_point =3D fe>get_xyz(); // The element shape functions evaluated at the quadrature points. const std::vector<std::vector<Real> >& phi =3D fe>get_phi(); // The element shape function gradients evaluated at the quadrature // points. const std::vector<std::vector<RealGradient> >& dphi =3D fe>get_dphi(); // A reference to the \p DofMap object for this system. The \p DofMap // object handles the index translation from node and element numbers // to degree of freedom numbers. const DofMap& dof_map =3D eigen_system.get_dof_map(); // The element mass and stiffness matrices. DenseMatrix<Number> Me; DenseMatrix<Number> Ke; // This vector will hold the degree of freedom indices for // the element. These define where in the global system // the element degrees of freedom get mapped. std::vector<unsigned int> dof_indices; // Now we will loop over all the elements in the mesh. // We will compute the element matrix contribution. MeshBase::const_element_iterator el =3D mesh.elements_begin(); const MeshBase::const_element_iterator end_el =3D mesh.elements_end(); for ( ; el !=3D end_el; ++el) { // Store a pointer to the element we are currently // working on. This allows for nicer syntax later. const Elem* elem =3D *el; // Get the degree of freedom indices for the // current element. These define where in the global // matrix and righthandside this element will // contribute to. dof_map.dof_indices (elem, dof_indices); // Compute the elementspecific data for the current // element. This involves computing the location of the // quadrature points (q_point) and the shape functions // (phi, dphi) for the current element. fe>reinit (elem); // Zero the element matrices before // summing them. We use the resize member here because // the number of degrees of freedom might have changed from // the last element. Note that this will be the case if the // element type is different (i.e. the last element was a // triangle, now we are on a quadrilateral). Ke.resize (dof_indices.size(), dof_indices.size()); Me.resize (dof_indices.size(), dof_indices.size()); // Now loop over the quadrature points. This handles // the numeric integration. // // We will build the element matrix. This involves // a double loop to integrate the test funcions (i) against // the trial functions (j). for (unsigned int qp=3D0; qp<qrule.n_points(); qp++) { const Real x =3D q_point[qp](0); =09for (unsigned int i=3D0; i<phi.size(); i++) =09 for (unsigned int j=3D0; j<phi.size(); j++) =09 { =09 Me(i,j) +=3D JxW[qp]*phi[i][qp]*phi[j][qp]; =09 Ke(i,j) +=3D .5*JxW[qp]*(dphi[i][qp]*dphi[j][qp] + pow(x, 2.)*phi[i][qp]*phi[j][qp]); =09 } } //Apply boundary conditions u(5) =3D u(5) =3D 0, follow the code in = ex0 double penalty =3D 1.e10; for(unsigned int s=3D0; s<elem>n_sides(); s++) { if(elem>neighbor(s) =3D=3D NULL) { AutoPtr<DofObject> node(elem>side(s)); for(unsigned int n=3D0; n<elem>n_nodes(); n++) { if(elem>node(n) =3D=3D node>id()) { Ke(n,n) +=3D penalty; Me(n,n) +=3D 0*penalty; } } } } // Finally, simply add the element contribution to the // overall matrices A and B. matrix_A.add_matrix (Ke, dof_indices); matrix_B.add_matrix (Me, dof_indices); } // end of element loop #endif // HAVE_SLEPC /** * All done! */ return; } =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D Thank you! David 
From: David Xu <dxu@my...>  20060304 17:59:32

Hi Everyone, I followed Steffen's new generalized eigen problem code in the new ex17 and wrote my first libmesh code to solve 1D timeindependent schrodinger equation with simple harmonic oscillator potential, the equation is 0.5u''+0.5*x^2*u=3DE*u . The analytical solution of the eigenvalue E are: 0.5, 1.5, 2.5, 3.5. 4.5 ... I set up the 1D domain [5 5] and apply the Dirichlet BC u(5)=3Du(5)=3D0 following ex0. The code compiled okay, but the eigenvalue results were huge numbers and I don't know what went wrong in my code. So, I'm attaching my code here, it should fairly straightforward since it's a simple 1D case. I'd really appreciate it if anyone can help me figure out the problem. Source code: =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D // 1D Time independent Schrodinger Equation Eigenvalue Solver // using simple harmonic potential // libMesh include files. #include "libmesh.h" #include "mesh.h" #include "mesh_generation.h" #include "gmv_io.h" #include "eigen_system.h" #include "equation_systems.h" #include "fe.h" #include "quadrature_gauss.h" #include "dense_matrix.h" #include "sparse_matrix.h" #include "numeric_vector.h" #include "dof_map.h" #include "edge_edge3.h" // Function prototype. This is the function that will assemble // the eigen system. Here, we will simply assemble a mass matrix. void assemble_mass(EquationSystems& es, const std::string& system_name); int main (int argc, char** argv) { // Initialize libMesh and the dependent libraries. libMesh::init (argc, argv); // This example is designed for the SLEPc eigen solver interface. #ifndef HAVE_SLEPC std::cerr << "ERROR: This example requires libMesh to be\n" << "compiled with SLEPc eigen solvers support!" << std::endl; return 0; #else // Braces are used to force object scope. { // Check for proper usage. if (argc < 3) { std::cerr << "\nUsage: " << argv[0] << " n <number of eigen values>" << std::endl; error(); } // Tell the user what we are doing. else { std::cout << "Running " << argv[0]; for (int i=3D1; i<argc; i++) std::cout << " " << argv[i]; std::cout << std::endl << std::endl; } // Set the dimensionality. const unsigned int dim =3D 1; // Get the number of eigen values to be computed from argv[2] const unsigned int nev =3D std::atoi(argv[2]); // Create a dimdimensional mesh. Mesh mesh (dim); // Use the internal mesh generator to create a uniform // grid on a 1D [5 5] domain. MeshTools::Generation::build_line (mesh, 100, 5., 5., EDGE3); // Print information about the mesh to the screen. mesh.print_info(); // Create an equation systems object. EquationSystems equation_systems (mesh); // Create a EigenSystem named "Eigensystem" and (for convenience) // use a reference to the system we create. EigenSystem & eigen_system =3D equation_systems.add_system<EigenSystem> ("Eigensystem"); // Declare the system variables. { // Adds the variable "p" to "Eigensystem". "p" // will be approximated using secondorder approximation. eigen_system.add_variable("u", SECOND); // Give the system a pointer to the matrix assembly // function defined below. eigen_system.attach_assemble_function (assemble_mass); // Set necessary parametrs used in EigenSystem::solve(), // i.e. the number of requested eigenpairs \p nev and the number // of basis vectors \p ncv used in the solution algorithm. Note that // ncv >=3D nev must hold and ncv >=3D 2*nev is recommended. equation_systems.parameters.set<unsigned int>("eigenpairs") =3D ne= v; equation_systems.parameters.set<unsigned int>("basis vectors") =3D ne= v*3; // Set the eigen solver type. SLEPc offers various solvers such as // the Arnoldi and subspace method. It // also offers interfaces to other solver packages (e.g. ARPACK). eigen_system.eigen_solver>set_eigensolver_type(ARNOLDI); // Set the solver tolerance and the maximum number of iterations. equation_systems.parameters.set<Real>("linear solver tolerance") =3D pow(TOLERANCE, 5./3.); equation_systems.parameters.set<unsigned int> ("linear solver maximum iterations") =3D 1000; // Set the type of the problem, here we deal with // a generalized Hermitian problem. eigen_system.set_eigenproblem_type(GHEP); // Set the eigenvalues to be computed. Note that not // all solvers support this. // eigen_system.eigen_solver>set_position_of_spectrum(SMALLEST_MAGNI= TUDE); // Initialize the data structures for the equation system. equation_systems.init(); // Prints information about the system to the screen. equation_systems.print_info(); } // Solve the system "Eigensystem". eigen_system.solve(); // Get the number of converged eigen pairs. unsigned int nconv =3D eigen_system.get_n_converged(); std::cout << "Number of converged eigenpairs: " << nconv << "\n" << std::endl; std::pair<Real, Real> myeig; // Get the last converged eigenpair if (nconv !=3D 0) { //print eigenvalues for(unsigned int n=3D0; n<=3Dnconv1; n++) { myeig=3Deigen_system.get_eigenpair(n); std::cout << "eig(" << n << ")=3D " << myeig.first << "\n" << std::= endl; } } else { std::cout << "WARNING: Solver did not converge!\n" << nconv << std::endl; } } #endif // HAVE_SLEPC // All done. return libMesh::close (); } void assemble_mass(EquationSystems& es, const std::string& system_name) { // It is a good idea to make sure we are assembling // the proper system. assert (system_name =3D=3D "Eigensystem"); #ifdef HAVE_SLEPC // Get a constant reference to the mesh object. const Mesh& mesh =3D es.get_mesh(); // The dimension that we are running. const unsigned int dim =3D mesh.mesh_dimension(); // Get a reference to our system. EigenSystem & eigen_system =3D es.get_system<EigenSystem> (system_name); // Get a constant reference to the Finite Element type // for the first (and only) variable in the system. FEType fe_type =3D eigen_system.get_dof_map().variable_type(0); // A reference to the two system matrices SparseMatrix<Number>& matrix_A =3D *eigen_system.matrix_A; SparseMatrix<Number>& matrix_B =3D *eigen_system.matrix_B; // Build a Finite Element object of the specified type. Since the // \p FEBase::build() member dynamically creates memory we will // store the object as an \p AutoPtr<FEBase>. This can be thought // of as a pointer that will clean up after itself. AutoPtr<FEBase> fe (FEBase::build(dim, fe_type)); // A Gauss quadrature rule for numerical integration. // Use the default quadrature order. //QGauss qrule (dim, fe_type.default_quadrature_order()); QGauss qrule(dim,FIFTH); // Tell the finite element object to use our quadrature rule. fe>attach_quadrature_rule (&qrule); // The element Jacobian * quadrature weight at each integration point. const std::vector<Real>& JxW =3D fe>get_JxW(); const std::vector<Point>& q_point =3D fe>get_xyz(); // The element shape functions evaluated at the quadrature points. const std::vector<std::vector<Real> >& phi =3D fe>get_phi(); // The element shape function gradients evaluated at the quadrature // points. const std::vector<std::vector<RealGradient> >& dphi =3D fe>get_dphi(); // A reference to the \p DofMap object for this system. The \p DofMap // object handles the index translation from node and element numbers // to degree of freedom numbers. const DofMap& dof_map =3D eigen_system.get_dof_map(); // The element mass and stiffness matrices. DenseMatrix<Number> Me; DenseMatrix<Number> Ke; // This vector will hold the degree of freedom indices for // the element. These define where in the global system // the element degrees of freedom get mapped. std::vector<unsigned int> dof_indices; // Now we will loop over all the elements in the mesh. // We will compute the element matrix contribution. MeshBase::const_element_iterator el =3D mesh.elements_begin(); const MeshBase::const_element_iterator end_el =3D mesh.elements_end(); for ( ; el !=3D end_el; ++el) { // Store a pointer to the element we are currently // working on. This allows for nicer syntax later. const Elem* elem =3D *el; // Get the degree of freedom indices for the // current element. These define where in the global // matrix and righthandside this element will // contribute to. dof_map.dof_indices (elem, dof_indices); // Compute the elementspecific data for the current // element. This involves computing the location of the // quadrature points (q_point) and the shape functions // (phi, dphi) for the current element. fe>reinit (elem); // Zero the element matrices before // summing them. We use the resize member here because // the number of degrees of freedom might have changed from // the last element. Note that this will be the case if the // element type is different (i.e. the last element was a // triangle, now we are on a quadrilateral). Ke.resize (dof_indices.size(), dof_indices.size()); Me.resize (dof_indices.size(), dof_indices.size()); // Now loop over the quadrature points. This handles // the numeric integration. // // We will build the element matrix. This involves // a double loop to integrate the test funcions (i) against // the trial functions (j). for (unsigned int qp=3D0; qp<qrule.n_points(); qp++) { const Real x =3D q_point[qp](0); for (unsigned int i=3D0; i<phi.size(); i++) for (unsigned int j=3D0; j<phi.size(); j++) { Me(i,j) +=3D JxW[qp]*phi[i][qp]*phi[j][qp]; Ke(i,j) +=3D .5*JxW[qp]*(dphi[i][qp]*dphi[j][qp] + pow(x, 2.)*phi[i][qp]*phi[j][qp]); } } //Apply boundary conditions u(5) =3D u(5) =3D 0, follow the code in = ex0 double penalty =3D 1.e10; for(unsigned int s=3D0; s<elem>n_sides(); s++) { if(elem>neighbor(s) =3D=3D NULL) { AutoPtr<DofObject> node(elem>side(s)); for(unsigned int n=3D0; n<elem>n_nodes(); n++) { if(elem>node(n) =3D=3D node>id()) { Ke(n,n) +=3D penalty; Me(n,n) +=3D 0*penalty; } } } } // Finally, simply add the element contribution to the // overall matrices A and B. matrix_A.add_matrix (Ke, dof_indices); matrix_B.add_matrix (Me, dof_indices); } // end of element loop #endif // HAVE_SLEPC /** * All done! */ return; } =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D Thank you! David 