25 #include <queso/Defines.h>
27 #ifdef QUESO_HAVE_LIBMESH_SLEPC
32 #include <queso/SharedPtr.h>
33 #include <queso/FunctionOperatorBuilder.h>
34 #include <queso/LibMeshNegativeLaplacianOperator.h>
35 #include <libmesh/libmesh_common.h>
36 #include <libmesh/mesh.h>
37 #include <libmesh/equation_systems.h>
38 #include <libmesh/condensed_eigen_system.h>
39 #include <libmesh/auto_ptr.h>
40 #include <libmesh/fe_base.h>
41 #include <libmesh/quadrature_gauss.h>
42 #include <libmesh/dense_matrix.h>
43 #include <libmesh/sparse_matrix.h>
44 #include <libmesh/dof_map.h>
45 #include <libmesh/id_types.h>
46 #include <libmesh/elem.h>
47 #include <libmesh/zero_function.h>
48 #include <libmesh/dirichlet_boundaries.h>
49 #include <libmesh/utility.h>
50 #include <libmesh/string_to_enum.h>
51 #include <libmesh/enum_order.h>
52 #include <libmesh/enum_fe_family.h>
56 LibMeshNegativeLaplacianOperator::LibMeshNegativeLaplacianOperator(
57 const FunctionOperatorBuilder & builder, libMesh::MeshBase & m)
58 : LibMeshOperatorBase(builder, m)
60 typename SharedPtr<libMesh::EquationSystems>::Type es(this->equation_systems);
64 libMesh::CondensedEigenSystem & eigen_system =
65 es->get_system<libMesh::CondensedEigenSystem>(
"Eigensystem");
70 unsigned int u_var = eigen_system.add_variable(
"u",
71 libMesh::Utility::string_to_enum<libMeshEnums::Order>(this->builder.order),
72 libMesh::Utility::string_to_enum<libMeshEnums::FEFamily>(this->builder.family));
76 eigen_system.attach_assemble_object(*
this);
80 es->parameters.set<
unsigned int>(
"eigenpairs") =
81 this->builder.num_req_eigenpairs;
82 es->parameters.set<
unsigned int>(
"basis vectors") =
83 this->builder.num_req_eigenpairs * 3;
86 es->parameters.set<libMesh::Real>(
"linear solver tolerance") =
87 pow(libMesh::TOLERANCE, 5./3.);
88 es->parameters.set<
unsigned int>(
"linear solver maximum iterations") = 1000;
92 eigen_system.set_eigenproblem_type(libMeshEnums::GHEP);
96 eigen_system.eigen_solver->set_position_of_spectrum
97 (libMeshEnums::SMALLEST_MAGNITUDE);
101 std::set<libMesh::boundary_id_type> boundary_ids;
102 boundary_ids.insert(0);
103 boundary_ids.insert(1);
104 boundary_ids.insert(2);
105 boundary_ids.insert(3);
108 std::vector<unsigned int> vars;
109 vars.push_back(u_var);
112 libMesh::ZeroFunction<> zf;
113 libMesh::DirichletBoundary dirichlet_bc(boundary_ids, vars, &zf);
114 eigen_system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);
120 std::set<unsigned int> global_dirichlet_dof_ids;
121 libMesh::DofConstraints::const_iterator
it =
122 eigen_system.get_dof_map().constraint_rows_begin();
124 for ( ; it != eigen_system.get_dof_map().constraint_rows_end(); ++
it) {
125 global_dirichlet_dof_ids.insert(it->first);
128 eigen_system.initialize_condensed_dofs(global_dirichlet_dof_ids);
131 eigen_system.solve();
134 this->nconv = eigen_system.get_n_converged();
137 LibMeshNegativeLaplacianOperator::~LibMeshNegativeLaplacianOperator()
141 void LibMeshNegativeLaplacianOperator::assemble()
143 #ifdef LIBMESH_HAVE_SLEPC
145 typename SharedPtr<libMesh::EquationSystems>::Type es(this->equation_systems);
148 const libMesh::MeshBase& mesh = es->get_mesh();
151 const unsigned int dim = mesh.mesh_dimension();
154 libMesh::EigenSystem & eigen_system = es->get_system<libMesh::EigenSystem>(
"Eigensystem");
158 libMesh::FEType fe_type = eigen_system.get_dof_map().variable_type(0);
161 libMesh::SparseMatrix<libMesh::Number>& matrix_A = *eigen_system.matrix_A;
162 libMesh::SparseMatrix<libMesh::Number>& matrix_B = *eigen_system.matrix_B;
168 libMesh::UniquePtr<libMesh::FEBase> fe (libMesh::FEBase::build(dim, fe_type));
172 libMesh::QGauss qrule(dim, fe_type.default_quadrature_order());
175 fe->attach_quadrature_rule (&qrule);
178 const std::vector<libMesh::Real>& JxW = fe->get_JxW();
181 const std::vector<std::vector<libMesh::Real> >& phi = fe->get_phi();
182 const std::vector<std::vector<libMesh::RealGradient> >& dphi = fe->get_dphi();
187 const libMesh::DofMap& dof_map = eigen_system.get_dof_map();
190 libMesh::DenseMatrix<libMesh::Number> Me;
191 libMesh::DenseMatrix<libMesh::Number> Ke;
196 std::vector<libMesh::dof_id_type> dof_indices;
204 libMesh::MeshBase::const_element_iterator el = mesh.active_local_elements_begin();
205 const libMesh::MeshBase::const_element_iterator end_el =
206 mesh.active_local_elements_end();
208 for ( ; el != end_el; ++el) {
211 const libMesh::Elem* elem = *el;
217 dof_map.dof_indices(elem, dof_indices);
231 Ke.resize(dof_indices.size(), dof_indices.size());
232 Me.resize(dof_indices.size(), dof_indices.size());
240 for (
unsigned int qp=0; qp < qrule.n_points(); qp++)
241 for (
unsigned int i=0; i < phi.size(); i++)
242 for (
unsigned int j=0; j < phi.size(); j++)
244 Me(i, j) += JxW[qp] * phi[i][qp] * phi[j][qp];
245 Ke(i, j) += JxW[qp] * (dphi[i][qp] * dphi[j][qp]);
262 matrix_A.add_matrix (Ke, dof_indices);
263 matrix_B.add_matrix (Me, dof_indices);
267 #endif // LIBMESH_HAVE_SLEPC
270 void LibMeshNegativeLaplacianOperator::print_info()
const
273 this->equation_systems->get_mesh().print_info();
274 this->equation_systems->print_info();
279 #endif // QUESO_HAVE_LIBMESH_SLEPC
that you receive source code or can get it if you want it