25 #include <queso/Defines.h>
27 #ifdef QUESO_HAVE_LIBMESH_SLEPC
32 #include <queso/FunctionOperatorBuilder.h>
33 #include <queso/LibMeshFunction.h>
34 #include <queso/LibMeshOperatorBase.h>
36 #include <libmesh/libmesh.h>
37 #include <libmesh/mesh.h>
38 #include <libmesh/mesh_generation.h>
39 #include <libmesh/equation_systems.h>
40 #include <libmesh/explicit_system.h>
41 #include <libmesh/eigen_system.h>
42 #include <libmesh/condensed_eigen_system.h>
43 #include <libmesh/numeric_vector.h>
44 #include <libmesh/exodusII_io.h>
45 #include <libmesh/parallel.h>
46 #include <libmesh/parallel_implementation.h>
50 LibMeshOperatorBase::LibMeshOperatorBase(
51 const FunctionOperatorBuilder & builder, libMesh::MeshBase & m)
53 equation_systems(new libMesh::EquationSystems(m)),
56 #ifndef LIBMESH_HAVE_SLEPC
57 if (m.processor_id() == 0)
58 std::cerr <<
"ERROR: This example requires libMesh to be\n"
59 <<
"compiled with SLEPc eigen solvers support!"
63 #ifdef LIBMESH_DEFAULT_SINGLE_PRECISION
65 libmesh_example_assert(
false,
"--disable-singleprecision");
68 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
71 libmesh_example_assert(
false,
"--disable-complex");
76 this->equation_systems->add_system<libMesh::CondensedEigenSystem>(
"Eigensystem");
78 #endif // LIBMESH_HAVE_SLEPC
81 LibMeshOperatorBase::~LibMeshOperatorBase()
85 void LibMeshOperatorBase::save_converged_evals(
const std::string & filename)
const
88 std::ofstream evals_file(filename.c_str());
90 std::pair<libMesh::Real, libMesh::Real> eval;
91 for (i = 0; i < this->nconv; i++) {
92 eval = this->equation_systems
93 ->get_system<libMesh::EigenSystem>(
"Eigensystem").get_eigenpair(i);
94 if (equation_systems->processor_id() == 0) {
95 evals_file << eval.first <<
" " << eval.second << std::endl;
98 if (this->equation_systems->processor_id() == 0) {
103 void LibMeshOperatorBase::save_converged_evec(
const std::string & filename,
104 unsigned int i)
const
106 if (i < this->nconv) {
107 typename SharedPtr<libMesh::EquationSystems>::Type es(this->equation_systems);
108 es->get_system<libMesh::EigenSystem>(
"Eigensystem").get_eigenpair(i);
109 libMesh::ExodusII_IO(es->get_mesh()).write_equation_systems(filename, *es);
112 std::cerr <<
"Warning: eigenpair " << i
113 <<
" did not converge. Not saving."
118 unsigned int LibMeshOperatorBase::get_num_converged()
const {
122 double LibMeshOperatorBase::get_eigenvalue(
unsigned int i)
const
124 if (i < this->nconv) {
125 std::pair<libMesh::Real, libMesh::Real> eval;
126 typename SharedPtr<libMesh::EquationSystems>::Type es(this->equation_systems);
127 eval = es->get_system<libMesh::EigenSystem>(
"Eigensystem").get_eigenpair(i);
135 double LibMeshOperatorBase::get_inverted_eigenvalue(
unsigned int i)
const
137 return 1.0 / this->get_eigenvalue(i);
140 libMesh::EquationSystems & LibMeshOperatorBase::get_equation_systems()
const
142 return *this->equation_systems;
145 typename SharedPtr<FunctionBase>::Type
146 LibMeshOperatorBase::inverse_kl_transform(std::vector<double> & xi,
150 typename SharedPtr<libMesh::EquationSystems>::Type es(this->equation_systems);
151 LibMeshFunction *kl =
new LibMeshFunction(this->builder, es->get_mesh());
156 this->equation_systems->comm().broadcast(xi);
158 typename SharedPtr<libMesh::EquationSystems>::Type kl_eq_sys(kl->get_equation_systems());
160 std::pair<libMesh::Real, libMesh::Real> eval;
161 for (i = 0; i < this->get_num_converged(); i++) {
162 eval = es->get_system<libMesh::EigenSystem>(
"Eigensystem").get_eigenpair(i);
163 kl_eq_sys->get_system<libMesh::ExplicitSystem>(
"Function").solution->add(
164 xi[i] / pow(eval.first, alpha / 2.0),
165 *es->get_system<libMesh::EigenSystem>(
"Eigensystem").solution);
168 typename SharedPtr<FunctionBase>::Type ap(kl);
174 #endif // QUESO_HAVE_LIBMESH_SLEPC