libMesh
optimization_ex2.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // <h1>Optimization Example 2 - Optimization with constraints</h1>
21 // \author David Knezevic
22 // \date 2015
23 //
24 // In this example we extend example 1 to demonstrate how to use
25 // OptimizationSystem's interface for imposing equality and inequality
26 // constraints.
27 
28 // C++ include files that we need
29 #include <iostream>
30 
31 // libMesh includes
32 #include "libmesh/libmesh.h"
33 #include "libmesh/mesh.h"
34 #include "libmesh/mesh_generation.h"
35 #include "libmesh/exodusII_io.h"
36 #include "libmesh/equation_systems.h"
37 #include "libmesh/fe.h"
38 #include "libmesh/quadrature_gauss.h"
39 #include "libmesh/dof_map.h"
40 #include "libmesh/sparse_matrix.h"
41 #include "libmesh/numeric_vector.h"
42 #include "libmesh/dense_matrix.h"
43 #include "libmesh/dense_vector.h"
44 #include "libmesh/elem.h"
45 #include "libmesh/boundary_info.h"
46 #include "libmesh/string_to_enum.h"
47 #include "libmesh/getpot.h"
48 #include "libmesh/zero_function.h"
49 #include "libmesh/dirichlet_boundaries.h"
50 #include "libmesh/optimization_system.h"
51 #include "libmesh/optimization_solver.h"
52 
53 // Bring in everything from the libMesh namespace
54 using namespace libMesh;
55 
69 {
70 private:
71 
75  OptimizationSystem & _sys;
76 
77 public:
78 
83 
89  void assemble_A_and_F();
90 
94  virtual Number objective (const NumericVector<Number> & soln,
95  OptimizationSystem & /*sys*/);
96 
100  virtual void gradient (const NumericVector<Number> & soln,
101  NumericVector<Number> & grad_f,
102  OptimizationSystem & /*sys*/);
103 
107  virtual void hessian (const NumericVector<Number> & soln,
108  SparseMatrix<Number> & H_f,
109  OptimizationSystem & /*sys*/);
110 
114  virtual void equality_constraints (const NumericVector<Number> & X,
115  NumericVector<Number> & C_eq,
116  OptimizationSystem & /*sys*/);
117 
121  virtual void equality_constraints_jacobian (const NumericVector<Number> & X,
122  SparseMatrix<Number> & C_eq_jac,
123  OptimizationSystem & /*sys*/);
124 
128  virtual void inequality_constraints (const NumericVector<Number> & X,
129  NumericVector<Number> & C_ineq,
130  OptimizationSystem & /*sys*/);
131 
135  virtual void inequality_constraints_jacobian (const NumericVector<Number> & X,
136  SparseMatrix<Number> & C_ineq_jac,
137  OptimizationSystem & /*sys*/);
138 
142  virtual void lower_and_upper_bounds (OptimizationSystem & sys);
143 
149  SparseMatrix<Number> * A_matrix;
150 
155  NumericVector<Number> * F_vector;
156 };
157 
159  _sys(sys_in)
160 {}
161 
163 {
164  A_matrix->zero();
165  F_vector->zero();
166 
167  const MeshBase & mesh = _sys.get_mesh();
168 
169  const unsigned int dim = mesh.mesh_dimension();
170  const unsigned int u_var = _sys.variable_number ("u");
171 
172  const DofMap & dof_map = _sys.get_dof_map();
173  FEType fe_type = dof_map.variable_type(u_var);
174  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
175  QGauss qrule (dim, fe_type.default_quadrature_order());
176  fe->attach_quadrature_rule (&qrule);
177 
178  const std::vector<Real> & JxW = fe->get_JxW();
179  const std::vector<std::vector<Real>> & phi = fe->get_phi();
180  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
181 
182  std::vector<dof_id_type> dof_indices;
183 
186 
187  for (const auto & elem : mesh.active_local_element_ptr_range())
188  {
189  dof_map.dof_indices (elem, dof_indices);
190 
191  const unsigned int n_dofs = dof_indices.size();
192 
193  fe->reinit (elem);
194 
195  Ke.resize (n_dofs, n_dofs);
196  Fe.resize (n_dofs);
197 
198  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
199  {
200  for (unsigned int dof_i=0; dof_i<n_dofs; dof_i++)
201  {
202  for (unsigned int dof_j=0; dof_j<n_dofs; dof_j++)
203  {
204  Ke(dof_i, dof_j) += JxW[qp] * (dphi[dof_j][qp]* dphi[dof_i][qp]);
205  }
206  Fe(dof_i) += JxW[qp] * phi[dof_i][qp];
207  }
208  }
209 
210  A_matrix->add_matrix (Ke, dof_indices);
211  F_vector->add_vector (Fe, dof_indices);
212  }
213 
214  A_matrix->close();
215  F_vector->close();
216 }
217 
219  OptimizationSystem & /*sys*/)
220 {
222 
223  A_matrix->vector_mult(*AxU, soln);
224  Number UTxAxU = AxU->dot(soln);
225 
226  Number UTxF = F_vector->dot(soln);
227 
228  return 0.5 * UTxAxU - UTxF;
229 }
230 
232  NumericVector<Number> & grad_f,
233  OptimizationSystem & /*sys*/)
234 {
235  grad_f.zero();
236 
237  A_matrix->vector_mult(grad_f, soln);
238  grad_f.add(-1, *F_vector);
239 }
240 
241 
243  SparseMatrix<Number> & H_f,
245 {
246  H_f.zero();
247  H_f.add(1., *A_matrix);
248 
249  // We also need to add the Hessian of the inequality and equality constraints,
250  // \sum_i^n_eq lambda_eq_i H_eq_i + \sum_i^n_ineq lambda_ineq_i H_ineq_i
251  // where lambda_eq and lambda_ineq denote Lagrange multipliers associated
252  // with the equality and inequality constraints, respectively.
253  // In this example, our equality constraints are linear, hence H_eq_i = 0.
254  // However, our inequality constraint is nonlinear, so it will contribute
255  // to the Hessian matrix.
256  sys.optimization_solver->get_dual_variables();
257  dof_id_type ineq_index = 0;
258  Number lambda_ineq_0 = 0.;
259  unsigned int lambda_rank = 0;
260  if ((sys.lambda_ineq->first_local_index() <= ineq_index) &&
261  (ineq_index < sys.lambda_ineq->last_local_index()))
262  {
263  lambda_ineq_0 = (*sys.lambda_ineq)(0);
264  lambda_rank = sys.comm().rank();
265  }
266 
267  // Sync lambda_rank across all processors.
268  sys.comm().sum(lambda_rank);
269  sys.comm().broadcast(lambda_rank, lambda_rank);
270 
271  if ((sys.get_dof_map().first_dof() <= 200) && (200 < sys.get_dof_map().end_dof()))
272  H_f.add(200, 200, 2. * lambda_ineq_0);
273 }
274 
276  NumericVector<Number> & C_eq,
277  OptimizationSystem & /*sys*/)
278 {
279  C_eq.zero();
280 
281  UniquePtr<NumericVector<Number>> X_localized =
283  X_localized->init(X.size(), false, SERIAL);
284  X.localize(*X_localized);
285 
286  std::vector<Number> constraint_values(3);
287  constraint_values[0] = (*X_localized)(17);
288  constraint_values[1] = (*X_localized)(23);
289  constraint_values[2] = (*X_localized)(98) + (*X_localized)(185);
290 
291  for (std::size_t i=0; i<constraint_values.size(); i++)
292  if ((C_eq.first_local_index() <= i) &&
293  (i < C_eq.last_local_index()))
294  C_eq.set(i, constraint_values[i]);
295 }
296 
298  SparseMatrix<Number> & C_eq_jac,
299  OptimizationSystem & sys)
300 {
301  C_eq_jac.zero();
302 
303  std::vector<std::vector<Number>> constraint_jac_values(3);
304  std::vector<std::vector<dof_id_type>> constraint_jac_indices(3);
305 
306  constraint_jac_values[0].resize(1);
307  constraint_jac_indices[0].resize(1);
308  constraint_jac_values[0][0] = 1.;
309  constraint_jac_indices[0][0] = 17;
310 
311  constraint_jac_values[1].resize(1);
312  constraint_jac_indices[1].resize(1);
313  constraint_jac_values[1][0] = 1.;
314  constraint_jac_indices[1][0] = 23;
315 
316  constraint_jac_values[2].resize(2);
317  constraint_jac_indices[2].resize(2);
318  constraint_jac_values[2][0] = 1.;
319  constraint_jac_values[2][1] = 1.;
320  constraint_jac_indices[2][0] = 98;
321  constraint_jac_indices[2][1] = 185;
322 
323  for (std::size_t i=0; i<constraint_jac_values.size(); i++)
324  for (std::size_t j=0; j<constraint_jac_values[i].size(); j++)
325  if ((sys.C_eq->first_local_index() <= i) &&
326  (i < sys.C_eq->last_local_index()))
327  {
328  dof_id_type col_index = constraint_jac_indices[i][j];
329  Number value = constraint_jac_values[i][j];
330  C_eq_jac.set(i, col_index, value);
331  }
332 }
333 
335  NumericVector<Number> & C_ineq,
336  OptimizationSystem & /*sys*/)
337 {
338  C_ineq.zero();
339 
340  UniquePtr<NumericVector<Number>> X_localized =
342  X_localized->init(X.size(), false, SERIAL);
343  X.localize(*X_localized);
344 
345  std::vector<Number> constraint_values(1);
346  constraint_values[0] = (*X_localized)(200)*(*X_localized)(200) + (*X_localized)(201) - 5.;
347 
348  for (std::size_t i=0; i<constraint_values.size(); i++)
349  if ((C_ineq.first_local_index() <= i) && (i < C_ineq.last_local_index()))
350  C_ineq.set(i, constraint_values[i]);
351 }
352 
354  SparseMatrix<Number> & C_ineq_jac,
355  OptimizationSystem & sys)
356 {
357  C_ineq_jac.zero();
358 
359  UniquePtr<NumericVector<Number>> X_localized =
361  X_localized->init(X.size(), false, SERIAL);
362  X.localize(*X_localized);
363 
364  std::vector<std::vector<Number>> constraint_jac_values(1);
365  std::vector<std::vector<dof_id_type>> constraint_jac_indices(1);
366 
367  constraint_jac_values[0].resize(2);
368  constraint_jac_indices[0].resize(2);
369  constraint_jac_values[0][0] = 2.* (*X_localized)(200);
370  constraint_jac_values[0][1] = 1.;
371  constraint_jac_indices[0][0] = 200;
372  constraint_jac_indices[0][1] = 201;
373 
374  for (std::size_t i=0; i<constraint_jac_values.size(); i++)
375  for (std::size_t j=0; j<constraint_jac_values[i].size(); j++)
376  if ((sys.C_ineq->first_local_index() <= i) &&
377  (i < sys.C_ineq->last_local_index()))
378  {
379  dof_id_type col_index = constraint_jac_indices[i][j];
380  Number value = constraint_jac_values[i][j];
381  C_ineq_jac.set(i, col_index, value);
382  }
383 
384 }
385 
387 {
388  for (unsigned int i=sys.get_dof_map().first_dof(); i<sys.get_dof_map().end_dof(); i++)
389  {
390  sys.get_vector("lower_bounds").set(i, -2.);
391  sys.get_vector("upper_bounds").set(i, 2.);
392  }
393 }
394 
395 int main (int argc, char ** argv)
396 {
397  LibMeshInit init (argc, argv);
398 
399 #ifndef LIBMESH_HAVE_PETSC_TAO
400 
401  libmesh_example_requires(false, "PETSc >= 3.5.0 with built-in TAO support");
402 
403 #elif LIBMESH_USE_COMPLEX_NUMBERS
404 
405  // According to
406  // http://www.mcs.anl.gov/research/projects/tao/documentation/installation.html
407  // TAO & PETSc-complex are currently mutually exclusive
408  libmesh_example_requires(false, "PETSc >= 3.5.0 with built-in TAO support & real-numbers only");
409 
410 #endif
411 
412  // We use a 2D domain.
413  libmesh_example_requires(LIBMESH_DIM > 1, "--disable-1D-only");
414 
415  // TAO is giving us problems in parallel?
416  if (init.comm().size() != 1)
417  {
418  libMesh::out << "This example can currently only be run in serial." << std::endl;
419  return 77;
420  }
421 
422  GetPot infile("optimization_ex2.in");
423  const std::string approx_order = infile("approx_order", "FIRST");
424  const std::string fe_family = infile("fe_family", "LAGRANGE");
425  const unsigned int n_elem = infile("n_elem", 10);
426 
427  Mesh mesh(init.comm());
429  n_elem,
430  n_elem,
431  -1., 1.,
432  -1., 1.,
433  QUAD9);
434 
435  mesh.print_info();
436 
437  EquationSystems equation_systems (mesh);
438 
439  OptimizationSystem & system =
440  equation_systems.add_system<OptimizationSystem> ("Optimization");
441 
442  // The default is to use PETSc/Tao solvers, but let the user change
443  // the optimization solver package on the fly.
444  {
445  const std::string optimization_solver_type = infile("optimization_solver_type",
446  "PETSC_SOLVERS");
447  SolverPackage sp = Utility::string_to_enum<SolverPackage>(optimization_solver_type);
450  system.optimization_solver.reset(new_solver.release());
451  }
452 
453  // Set tolerances and maximum iteration counts directly on the optimization solver.
454  system.optimization_solver->max_objective_function_evaluations = 128;
455  system.optimization_solver->objective_function_relative_tolerance = 1.e-4;
456  system.optimization_solver->verbose = true;
457 
458  AssembleOptimization assemble_opt(system);
459 
460  system.add_variable("u",
461  Utility::string_to_enum<Order> (approx_order),
462  Utility::string_to_enum<FEFamily>(fe_family));
463 
464  system.optimization_solver->objective_object = &assemble_opt;
465  system.optimization_solver->gradient_object = &assemble_opt;
466  system.optimization_solver->hessian_object = &assemble_opt;
467  system.optimization_solver->equality_constraints_object = &assemble_opt;
468  system.optimization_solver->equality_constraints_jacobian_object = &assemble_opt;
469  system.optimization_solver->inequality_constraints_object = &assemble_opt;
470  system.optimization_solver->inequality_constraints_jacobian_object = &assemble_opt;
471  system.optimization_solver->lower_and_upper_bounds_object = &assemble_opt;
472 
473  // system.matrix and system.rhs are used for the gradient and Hessian,
474  // so in this case we add an extra matrix and vector to store A and F.
475  // This makes it easy to write the code for evaluating the objective,
476  // gradient, and hessian.
477  system.add_matrix("A_matrix");
478  system.add_vector("F_vector");
479  assemble_opt.A_matrix = &system.get_matrix("A_matrix");
480  assemble_opt.F_vector = &system.get_vector("F_vector");
481 
482 
483  equation_systems.init();
484  equation_systems.print_info();
485 
486  assemble_opt.assemble_A_and_F();
487 
488  {
489  std::vector<std::set<numeric_index_type>> constraint_jac_sparsity;
490  std::set<numeric_index_type> sparsity_row;
491  sparsity_row.insert(17);
492  constraint_jac_sparsity.push_back(sparsity_row);
493  sparsity_row.clear();
494 
495  sparsity_row.insert(23);
496  constraint_jac_sparsity.push_back(sparsity_row);
497  sparsity_row.clear();
498 
499  sparsity_row.insert(98);
500  sparsity_row.insert(185);
501  constraint_jac_sparsity.push_back(sparsity_row);
502 
503  system.initialize_equality_constraints_storage(constraint_jac_sparsity);
504  }
505 
506  {
507  std::vector<std::set<numeric_index_type>> constraint_jac_sparsity;
508  std::set<numeric_index_type> sparsity_row;
509  sparsity_row.insert(200);
510  sparsity_row.insert(201);
511  constraint_jac_sparsity.push_back(sparsity_row);
512 
513  system.initialize_inequality_constraints_storage(constraint_jac_sparsity);
514  }
515 
516  // We need to close the matrix so that we can use it to store the
517  // Hessian during the solve.
518  system.matrix->close();
519  system.solve();
520 
521  // Print convergence information
522  system.optimization_solver->print_converged_reason();
523 
524  // Write the solution to file if the optimization solver converged
525  if (system.optimization_solver->get_converged_reason() > 0)
526  {
527  std::stringstream filename;
528  ExodusII_IO (mesh).write_equation_systems("optimization_soln.exo",
529  equation_systems);
530  }
531 
532  return 0;
533 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
Abstract base class to be used to calculate the inequality constraints.
OptimizationSystem & _sys
Keep a reference to the OptimizationSystem.
This is the EquationSystems class.
virtual void lower_and_upper_bounds(OptimizationSystem &sys)
Evaluate the lower and upper bounds vectors.
virtual Number objective(const NumericVector< Number > &soln, OptimizationSystem &)
Evaluate the objective function.
void initialize_inequality_constraints_storage(const std::vector< std::set< numeric_index_type >> &constraint_jac_sparsity)
Initialize storage for the inequality constraints, as per initialize_equality_constraints_storage.
virtual numeric_index_type size() const =0
virtual numeric_index_type last_local_index() const =0
virtual void init(const numeric_index_type, const numeric_index_type, const bool=false, const ParallelType=AUTOMATIC)=0
Change the dimension of the vector to N.
Abstract base class to be used to calculate the objective function for optimization.
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:1697
unsigned int size() const
Definition: parallel.h:726
dof_id_type n_elem(const MeshBase::const_element_iterator &begin, const MeshBase::const_element_iterator &end)
Count up the number of elements of a specific type (as defined by an iterator range).
Definition: mesh_tools.C:656
AssembleOptimization(OptimizationSystem &sys_in)
Constructor.
UniquePtr< NumericVector< Number > > lambda_ineq
unsigned int dim
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
Definition: exodusII_io.h:52
ImplicitSystem & sys
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:350
virtual void add_vector(const T *v, const std::vector< numeric_index_type > &dof_indices)
Computes , where v is a pointer and each dof_indices[i] specifies where to add value v[i]...
unsigned int add_variable(const std::string &var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=libmesh_nullptr)
Adds the variable var to the list of variables for this system.
Definition: system.C:1101
SolverPackage
Defines an enum for various linear solver packages.
virtual void solve() libmesh_override
Solves the optimization problem.
SparseMatrix< Number > & add_matrix(const std::string &mat_name)
Adds the additional matrix mat_name to this system.
MeshBase & mesh
UniquePtr< OptimizationSolver< Number > > optimization_solver
The OptimizationSolver that is used for performing the optimization.
NumericVector< Number > * F_vector
Vector for storing F.
void vector_mult(NumericVector< T > &dest, const NumericVector< T > &arg) const
Multiplies the matrix by the NumericVector arg and stores the result in NumericVector dest...
virtual void equality_constraints(const NumericVector< Number > &X, NumericVector< Number > &C_eq, OptimizationSystem &)
Evaluate the equality constraints.
void build_square(UnstructuredMesh &mesh, const unsigned int nx, const unsigned int ny, const Real xmin=0., const Real xmax=1., const Real ymin=0., const Real ymax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
A specialized build_cube() for 2D meshes.
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
Definition: libmesh.h:62
The libMesh namespace provides an interface to certain functionality in the library.
virtual void equality_constraints_jacobian(const NumericVector< Number > &X, SparseMatrix< Number > &C_eq_jac, OptimizationSystem &)
Evaluate the equality constraints Jacobian.
This is the MeshBase class.
Definition: mesh_base.h:68
virtual void zero()=0
Set all entries to zero.
virtual void set(const numeric_index_type i, const numeric_index_type j, const T value)=0
Set the element (i,j) to value.
virtual void add(const numeric_index_type i, const numeric_index_type j, const T value)=0
Add value to the element (i,j).
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
virtual SimpleRange< element_iterator > active_local_element_ptr_range()=0
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
dof_id_type first_dof(const processor_id_type proc) const
Definition: dof_map.h:535
SparseMatrix< Number > * A_matrix
Sparse matrix for storing the matrix A.
unsigned short int variable_number(const std::string &var) const
Definition: system.C:1263
dof_id_type end_dof(const processor_id_type proc) const
Definition: dof_map.h:577
virtual void add_matrix(const DenseMatrix< T > &dm, const std::vector< numeric_index_type > &rows, const std::vector< numeric_index_type > &cols)=0
Add the full matrix dm to the SparseMatrix.
const MeshBase & get_mesh() const
Definition: system.h:2014
This class encapsulate all functionality required for assembling the objective function, gradient, and hessian.
NumericVector< Number > & add_vector(const std::string &vec_name, const bool projections=true, const ParallelType type=PARALLEL)
Adds the additional vector vec_name to this system.
Definition: system.C:681
const DofMap & get_dof_map() const
Definition: system.h:2030
virtual void zero()=0
Set all entries to 0.
Abstract base class to be used to calculate the equality constraints.
virtual void inequality_constraints_jacobian(const NumericVector< Number > &X, SparseMatrix< Number > &C_ineq_jac, OptimizationSystem &)
Evaluate the inequality constraints Jacobian.
void print_info(std::ostream &os=libMesh::out) const
Prints information about the equation systems, by default to libMesh::out.
Order default_quadrature_order() const
Definition: fe_type.h:332
This System subclass enables us to assemble an objective function, gradient, Hessian and bounds for o...
int main(int argc, char **argv)
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
virtual System & add_system(const std::string &system_type, const std::string &name)
Add the system of type system_type named name to the systems array.
const Parallel::Communicator & comm() const
Definition: libmesh.h:81
void broadcast(T &data, const unsigned int root_id=0) const
Take a local value and broadcast it to all processors.
const NumericVector< Number > & get_vector(const std::string &vec_name) const
Definition: system.C:794
virtual void close()=0
Calls the NumericVector&#39;s internal assembly routines, ensuring that the values are consistent across ...
const SparseMatrix< Number > & get_matrix(const std::string &mat_name) const
Abstract base class to be used to calculate the gradient of an objective function.
virtual void hessian(const NumericVector< Number > &soln, SparseMatrix< Number > &H_f, OptimizationSystem &)
Evaluate the Hessian.
virtual void write_equation_systems(const std::string &, const EquationSystems &, const std::set< std::string > *system_names=libmesh_nullptr)
This method implements writing a mesh with data to a specified file where the data is taken from the ...
Definition: mesh_output.C:31
Abstract base class to be used to calculate the Jacobian of the equality constraints.
virtual void close()=0
Calls the SparseMatrix&#39;s internal assembly routines, ensuring that the values are consistent across p...
void initialize_equality_constraints_storage(const std::vector< std::set< numeric_index_type >> &constraint_jac_sparsity)
Initialize storage for the equality constraints, and the corresponding Jacobian.
This base class can be inherited from to provide interfaces to optimization solvers from different pa...
virtual UniquePtr< NumericVector< T > > zero_clone() const =0
const Parallel::Communicator & comm() const
OStreamProxy out
SparseMatrix< Number > * matrix
The system matrix.
virtual void gradient(const NumericVector< Number > &soln, NumericVector< Number > &grad_f, OptimizationSystem &)
Evaluate the gradient.
unsigned int mesh_dimension() const
Definition: mesh_base.C:148
Abstract base class to be used to calculate the Hessian of an objective function. ...
static const bool value
Definition: xdr_io.C:108
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:776
This class implements specific orders of Gauss quadrature.
virtual numeric_index_type first_local_index() const =0
virtual T dot(const NumericVector< T > &v) const =0
unsigned int rank() const
Definition: parallel.h:724
virtual void set(const numeric_index_type i, const T value)=0
Sets v(i) = value.
virtual void init()
Initialize all the systems.
virtual void inequality_constraints(const NumericVector< Number > &X, NumericVector< Number > &C_ineq, OptimizationSystem &)
Evaluate the inequality constraints.
void assemble_A_and_F()
The optimization problem we consider here is: min_U 0.5*U^T A U - U^T F.
virtual void add(const numeric_index_type i, const T value)=0
Adds value to each entry of the vector.
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
UniquePtr< NumericVector< Number > > C_eq
The vector that stores equality constraints.
void sum(T &r) const
Take a local variable and replace it with the sum of it&#39;s values on all processors.
Abstract base class to be used to calculate the lower and upper bounds for all dofs in the system...
void print_info(std::ostream &os=libMesh::out) const
Prints relevant information about the mesh.
Definition: mesh_base.C:448
virtual void localize(std::vector< T > &v_local) const =0
Creates a copy of the global vector in the local vector v_local.
uint8_t dof_id_type
Definition: id_types.h:64
UniquePtr< NumericVector< Number > > C_ineq
The vector that stores inequality constraints.
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Fills the vector di with the global degree of freedom indices for the element.
Definition: dof_map.C:1917
Abstract base class to be used to calculate the Jacobian of the inequality constraints.