25 #include <queso/Defines.h>
27 #ifdef QUESO_HAVE_LIBMESH
32 #include <boost/shared_ptr.hpp>
34 #include <queso/FunctionOperatorBuilder.h>
35 #include <queso/LibMeshNegativeLaplacianOperator.h>
36 #include <libmesh/libmesh_common.h>
37 #include <libmesh/mesh.h>
38 #include <libmesh/equation_systems.h>
39 #include <libmesh/condensed_eigen_system.h>
40 #include <libmesh/auto_ptr.h>
41 #include <libmesh/fe_base.h>
42 #include <libmesh/quadrature_gauss.h>
43 #include <libmesh/dense_matrix.h>
44 #include <libmesh/sparse_matrix.h>
45 #include <libmesh/dof_map.h>
46 #include <libmesh/id_types.h>
47 #include <libmesh/elem.h>
48 #include <libmesh/zero_function.h>
49 #include <libmesh/dirichlet_boundaries.h>
50 #include <libmesh/utility.h>
51 #include <libmesh/string_to_enum.h>
52 #include <libmesh/enum_order.h>
53 #include <libmesh/enum_fe_family.h>
57 LibMeshNegativeLaplacianOperator::LibMeshNegativeLaplacianOperator(
58 const FunctionOperatorBuilder & builder, libMesh::MeshBase & m)
59 : LibMeshOperatorBase(builder, m)
61 boost::shared_ptr<libMesh::EquationSystems> es(this->equation_systems);
65 libMesh::CondensedEigenSystem & eigen_system =
66 es->get_system<libMesh::CondensedEigenSystem>(
"Eigensystem");
71 unsigned int u_var = eigen_system.add_variable(
"u",
72 libMesh::Utility::string_to_enum<libMeshEnums::Order>(this->builder.order),
73 libMesh::Utility::string_to_enum<libMeshEnums::FEFamily>(this->builder.family));
77 eigen_system.attach_assemble_object(*
this);
81 es->parameters.set<
unsigned int>(
"eigenpairs") =
82 this->builder.num_req_eigenpairs;
83 es->parameters.set<
unsigned int>(
"basis vectors") =
84 this->builder.num_req_eigenpairs * 3;
87 es->parameters.set<libMesh::Real>(
"linear solver tolerance") = pow(TOLERANCE, 5./3.);
88 es->parameters.set<
unsigned int>(
"linear solver maximum iterations") = 1000;
92 eigen_system.set_eigenproblem_type(GHEP);
96 eigen_system.eigen_solver->set_position_of_spectrum(SMALLEST_MAGNITUDE);
100 std::set<libMesh::boundary_id_type> boundary_ids;
101 boundary_ids.insert(0);
102 boundary_ids.insert(1);
103 boundary_ids.insert(2);
104 boundary_ids.insert(3);
107 std::vector<unsigned int> vars;
108 vars.push_back(u_var);
111 libMesh::ZeroFunction<> zf;
112 libMesh::DirichletBoundary dirichlet_bc(boundary_ids, vars, &zf);
113 eigen_system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);
119 std::set<unsigned int> global_dirichlet_dof_ids;
120 libMesh::DofConstraints::const_iterator it =
121 eigen_system.get_dof_map().constraint_rows_begin();
123 for ( ; it != eigen_system.get_dof_map().constraint_rows_end(); ++it) {
124 global_dirichlet_dof_ids.insert(it->first);
127 eigen_system.initialize_condensed_dofs(global_dirichlet_dof_ids);
130 eigen_system.solve();
133 this->nconv = eigen_system.get_n_converged();
136 LibMeshNegativeLaplacianOperator::~LibMeshNegativeLaplacianOperator()
140 void LibMeshNegativeLaplacianOperator::assemble()
142 #ifdef LIBMESH_HAVE_SLEPC
144 boost::shared_ptr<libMesh::EquationSystems> es(this->equation_systems);
147 const libMesh::MeshBase& mesh = es->get_mesh();
150 const unsigned int dim = mesh.mesh_dimension();
153 libMesh::EigenSystem & eigen_system = es->get_system<libMesh::EigenSystem>(
"Eigensystem");
157 libMesh::FEType fe_type = eigen_system.get_dof_map().variable_type(0);
160 libMesh::SparseMatrix<libMesh::Number>& matrix_A = *eigen_system.matrix_A;
161 libMesh::SparseMatrix<libMesh::Number>& matrix_B = *eigen_system.matrix_B;
167 libMesh::AutoPtr<libMesh::FEBase> fe (libMesh::FEBase::build(dim, fe_type));
171 libMesh::QGauss qrule(dim, fe_type.default_quadrature_order());
174 fe->attach_quadrature_rule (&qrule);
177 const std::vector<libMesh::Real>& JxW = fe->get_JxW();
180 const std::vector<std::vector<libMesh::Real> >& phi = fe->get_phi();
181 const std::vector<std::vector<libMesh::RealGradient> >& dphi = fe->get_dphi();
186 const libMesh::DofMap& dof_map = eigen_system.get_dof_map();
189 libMesh::DenseMatrix<libMesh::Number> Me;
190 libMesh::DenseMatrix<libMesh::Number> Ke;
195 std::vector<libMesh::dof_id_type> dof_indices;
203 libMesh::MeshBase::const_element_iterator el = mesh.active_local_elements_begin();
204 const libMesh::MeshBase::const_element_iterator end_el =
205 mesh.active_local_elements_end();
207 for ( ; el != end_el; ++el) {
210 const libMesh::Elem* elem = *el;
216 dof_map.dof_indices(elem, dof_indices);
230 Ke.resize(dof_indices.size(), dof_indices.size());
231 Me.resize(dof_indices.size(), dof_indices.size());
239 for (
unsigned int qp=0; qp < qrule.n_points(); qp++)
240 for (
unsigned int i=0; i < phi.size(); i++)
241 for (
unsigned int j=0; j < phi.size(); j++)
243 Me(i, j) += JxW[qp] * phi[i][qp] * phi[j][qp];
244 Ke(i, j) += JxW[qp] * (dphi[i][qp] * dphi[j][qp]);
261 matrix_A.add_matrix (Ke, dof_indices);
262 matrix_B.add_matrix (Me, dof_indices);
266 #endif // LIBMESH_HAVE_SLEPC
269 void LibMeshNegativeLaplacianOperator::print_info()
const
272 this->equation_systems->get_mesh().print_info();
273 this->equation_systems->print_info();
278 #endif // QUESO_HAVE_LIBMESH