queso-0.57.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
LibMeshNegativeLaplacianOperator.C
Go to the documentation of this file.
1 //-----------------------------------------------------------------------bl-
2 //--------------------------------------------------------------------------
3 //
4 // QUESO - a library to support the Quantification of Uncertainty
5 // for Estimation, Simulation and Optimization
6 //
7 // Copyright (C) 2008-2017 The PECOS Development Team
8 //
9 // This library is free software; you can redistribute it and/or
10 // modify it under the terms of the Version 2.1 GNU Lesser General
11 // Public License as published by the Free Software Foundation.
12 //
13 // This library is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // Lesser General Public License for more details.
17 //
18 // You should have received a copy of the GNU Lesser General Public
19 // License along with this library; if not, write to the Free Software
20 // Foundation, Inc. 51 Franklin Street, Fifth Floor,
21 // Boston, MA 02110-1301 USA
22 //
23 //-----------------------------------------------------------------------el-
24 
25 #include <queso/Defines.h>
26 
27 #ifdef QUESO_HAVE_LIBMESH_SLEPC
28 
29 #include <set>
30 #include <vector>
31 
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>
53 
54 namespace QUESO {
55 
57  const FunctionOperatorBuilder & builder, libMesh::MeshBase & m)
58  : LibMeshOperatorBase(builder, m)
59 {
61 
62  // Give the system a pointer to the matrix assembly
63  // function defined below.
64  libMesh::CondensedEigenSystem & eigen_system =
65  es->get_system<libMesh::CondensedEigenSystem>("Eigensystem");
66 
67  // Declare the system variables.
68  // Adds the variable "u" to "Eigensystem". "u"
69  // will be approximated using second-order approximation.
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));
73 
74  // This works because *this is a subclass of System::Assembly
75  // and requires the class to implement an 'assemble'
76  eigen_system.attach_assemble_object(*this);
77 
78  // Set the number of requested eigenpairs \p n_evals and the number
79  // of basis vectors used in the solution algorithm.
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;
84 
85  // Set the solver tolerance and the maximum number of iterations.
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;
89 
90  // Set the type of the problem, here we deal with
91  // a generalized Hermitian problem.
92  eigen_system.set_eigenproblem_type(libMeshEnums::GHEP);
93 
94  // Order the eigenvalues "smallest first"
95  // This hoses performance?
96  eigen_system.eigen_solver->set_position_of_spectrum
97  (libMeshEnums::SMALLEST_MAGNITUDE);
98 
99  // Set up the boundary (only works if this->m is a square)
100  // We'll just the whole boundary to be Dirichlet, because why not
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);
106 
107  // Assign which variables must satisfy the boundary constraints
108  std::vector<unsigned int> vars;
109  vars.push_back(u_var);
110 
111  // Add the boundary object to the eigensystem
112  libMesh::ZeroFunction<> zf;
113  libMesh::DirichletBoundary dirichlet_bc(boundary_ids, vars, &zf);
114  eigen_system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);
115 
116  // Initialize the data structures for the equation system.
117  es->init();
118 
119  // Here we obtain the ids of the Dirichlet dofs
120  std::set<unsigned int> global_dirichlet_dof_ids;
121  libMesh::DofConstraints::const_iterator it =
122  eigen_system.get_dof_map().constraint_rows_begin();
123 
124  for ( ; it != eigen_system.get_dof_map().constraint_rows_end(); ++it) {
125  global_dirichlet_dof_ids.insert(it->first);
126  }
127 
128  eigen_system.initialize_condensed_dofs(global_dirichlet_dof_ids);
129 
130  // Solve the system "Eigensystem".
131  eigen_system.solve();
132 
133  // Get the number of converged eigen pairs.
134  this->nconv = eigen_system.get_n_converged();
135 }
136 
138 {
139 }
140 
142 {
143 #ifdef LIBMESH_HAVE_SLEPC
144 
146 
147  // Get a constant reference to the mesh object.
148  const libMesh::MeshBase& mesh = es->get_mesh();
149 
150  // The dimension that we are running.
151  const unsigned int dim = mesh.mesh_dimension();
152 
153  // Get a reference to our system.
154  libMesh::EigenSystem & eigen_system = es->get_system<libMesh::EigenSystem>("Eigensystem");
155 
156  // Get a constant reference to the Finite Element type
157  // for the first (and only) variable in the system.
158  libMesh::FEType fe_type = eigen_system.get_dof_map().variable_type(0);
159 
160  // A reference to the two system matrices
161  libMesh::SparseMatrix<libMesh::Number>& matrix_A = *eigen_system.matrix_A;
162  libMesh::SparseMatrix<libMesh::Number>& matrix_B = *eigen_system.matrix_B;
163 
164  // Build a Finite Element object of the specified type. Since the
165  // \p FEBase::build() member dynamically creates memory we will
166  // store the object as an \p UniquePtr<FEBase>. This can be thought
167  // of as a pointer that will clean up after itself.
168  libMesh::UniquePtr<libMesh::FEBase> fe (libMesh::FEBase::build(dim, fe_type));
169 
170  // A Gauss quadrature rule for numerical integration.
171  // Use the default quadrature order.
172  libMesh::QGauss qrule(dim, fe_type.default_quadrature_order());
173 
174  // Tell the finite element object to use our quadrature rule.
175  fe->attach_quadrature_rule (&qrule);
176 
177  // The element Jacobian * quadrature weight at each integration point.
178  const std::vector<libMesh::Real>& JxW = fe->get_JxW();
179 
180  // The element shape functions evaluated at the quadrature points.
181  const std::vector<std::vector<libMesh::Real> >& phi = fe->get_phi();
182  const std::vector<std::vector<libMesh::RealGradient> >& dphi = fe->get_dphi();
183 
184  // A reference to the \p DofMap object for this system. The \p DofMap
185  // object handles the index translation from node and element numbers
186  // to degree of freedom numbers.
187  const libMesh::DofMap& dof_map = eigen_system.get_dof_map();
188 
189  // The element stiffness matrix.
190  libMesh::DenseMatrix<libMesh::Number> Me;
191  libMesh::DenseMatrix<libMesh::Number> Ke;
192 
193  // This vector will hold the degree of freedom indices for
194  // the element. These define where in the global system
195  // the element degrees of freedom get mapped.
196  std::vector<libMesh::dof_id_type> dof_indices;
197 
198  // Now we will loop over all the elements in the mesh that
199  // live on the local processor. We will compute the element
200  // matrix and right-hand-side contribution. In case users
201  // later modify this program to include refinement, we will
202  // be safe and will only consider the active elements;
203  // hence we use a variant of the \p active_elem_iterator.
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();
207 
208  for ( ; el != end_el; ++el) {
209  // Store a pointer to the element we are currently
210  // working on. This allows for nicer syntax later.
211  const libMesh::Elem* elem = *el;
212 
213  // Get the degree of freedom indices for the
214  // current element. These define where in the global
215  // matrix and right-hand-side this element will
216  // contribute to.
217  dof_map.dof_indices(elem, dof_indices);
218 
219  // Compute the element-specific data for the current
220  // element. This involves computing the location of the
221  // quadrature points (q_point) and the shape functions
222  // (phi, dphi) for the current element.
223  fe->reinit(elem);
224 
225  // Zero the element matrices before
226  // summing them. We use the resize member here because
227  // the number of degrees of freedom might have changed from
228  // the last element. Note that this will be the case if the
229  // element type is different (i.e. the last element was a
230  // triangle, now we are on a quadrilateral).
231  Ke.resize(dof_indices.size(), dof_indices.size());
232  Me.resize(dof_indices.size(), dof_indices.size());
233 
234  // Now loop over the quadrature points. This handles
235  // the numeric integration.
236  //
237  // We will build the element matrix. This involves
238  // a double loop to integrate the test funcions (i) against
239  // the trial functions (j).
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++)
243  {
244  Me(i, j) += JxW[qp] * phi[i][qp] * phi[j][qp];
245  Ke(i, j) += JxW[qp] * (dphi[i][qp] * dphi[j][qp]);
246  }
247 
248  // On an unrefined mesh, constrain_element_matrix does
249  // nothing. If this assembly function is ever repurposed to
250  // run on a refined mesh, getting the hanging node constraints
251  // right will be important. Note that, even with
252  // asymmetric_constraint_rows = false, the constrained dof
253  // diagonals still exist in the matrix, with diagonal entries
254  // that are there to ensure non-singular matrices for linear
255  // solves but which would generate positive non-physical
256  // eigenvalues for eigensolves.
257  // dof_map.constrain_element_matrix(Ke, dof_indices, false);
258  // dof_map.constrain_element_matrix(Me, dof_indices, false);
259 
260  // Finally, simply add the element contribution to the
261  // overall matrices A and B.
262  matrix_A.add_matrix (Ke, dof_indices);
263  matrix_B.add_matrix (Me, dof_indices);
264 
265  } // end of element loop
266 
267 #endif // LIBMESH_HAVE_SLEPC
268 }
269 
271 {
272  // Prints information about the system to the screen.
273  this->equation_systems->get_mesh().print_info();
274  this->equation_systems->print_info();
275 }
276 
277 } // End namespace QUESO
278 
279 #endif // QUESO_HAVE_LIBMESH_SLEPC
std::string order
String to store the polynomial order to use. Default is &quot;FIRST&quot;.
virtual void assemble()
Method to assemble the mass and stiffness matrices.
Abstract base class for operator objects using libmesh in the backend.
LibMeshNegativeLaplacianOperator(const FunctionOperatorBuilder &builder, libMesh::MeshBase &m)
Construct the negative Laplacian operator on the libmesh mesh m with builder builder.
virtual void print_info() const
Print libmesh related foo.
SharedPtr< libMesh::EquationSystems >::Type equation_systems
unsigned int nconv
The number of converged eigenpairs.
int dim
Definition: ann_test.cpp:472
unsigned int num_req_eigenpairs
Number of eigenpairs to request when building an operator. Default is 0.

Generated on Tue Jun 5 2018 19:48:55 for queso-0.57.1 by  doxygen 1.8.5