MFEM  v4.3.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
ex19.cpp
Go to the documentation of this file.
1 // MFEM Example 19
2 //
3 // Compile with: make ex19
4 //
5 // Sample runs:
6 // ex19 -m ../data/beam-quad.mesh
7 // ex19 -m ../data/beam-tri.mesh
8 // ex19 -m ../data/beam-hex.mesh
9 // ex19 -m ../data/beam-tet.mesh
10 // ex19 -m ../data/beam-wedge.mesh
11 // ex19 -m ../data/beam-quad-amr.mesh
12 //
13 // Description: This examples solves a quasi-static incompressible nonlinear
14 // elasticity problem of the form 0 = H(x), where H is an
15 // incompressible hyperelastic model and x is a block state vector
16 // containing displacement and pressure variables. The geometry of
17 // the domain is assumed to be as follows:
18 //
19 // +---------------------+
20 // boundary --->| |<--- boundary
21 // attribute 1 | | attribute 2
22 // (fixed) +---------------------+ (fixed, nonzero)
23 //
24 // The example demonstrates the use of block nonlinear operators
25 // (the class RubberOperator defining H(x)) as well as a nonlinear
26 // Newton solver for the quasi-static problem. Each Newton step
27 // requires the inversion of a Jacobian matrix, which is done
28 // through a (preconditioned) inner solver. The specialized block
29 // preconditioner is implemented as a user-defined solver.
30 //
31 // We recommend viewing examples 2, 5, and 10 before viewing this
32 // example.
33 
34 #include "mfem.hpp"
35 #include <memory>
36 #include <iostream>
37 #include <fstream>
38 
39 using namespace std;
40 using namespace mfem;
41 
42 class GeneralResidualMonitor : public IterativeSolverMonitor
43 {
44 public:
45  GeneralResidualMonitor(const std::string& prefix_, int print_lvl)
46  : prefix(prefix_)
47  {
48  print_level = print_lvl;
49  }
50 
51  virtual void MonitorResidual(int it, double norm, const Vector &r, bool final);
52 
53 private:
54  const std::string prefix;
55  int print_level;
56  mutable double norm0;
57 };
58 
59 void GeneralResidualMonitor::MonitorResidual(int it, double norm,
60  const Vector &r, bool final)
61 {
62  if (print_level == 1 || (print_level == 3 && (final || it == 0)))
63  {
64  mfem::out << prefix << " iteration " << setw(2) << it
65  << " : ||r|| = " << norm;
66  if (it > 0)
67  {
68  mfem::out << ", ||r||/||r_0|| = " << norm/norm0;
69  }
70  else
71  {
72  norm0 = norm;
73  }
74  mfem::out << '\n';
75  }
76 }
77 
78 // Custom block preconditioner for the Jacobian of the incompressible nonlinear
79 // elasticity operator. It has the form
80 //
81 // P^-1 = [ K^-1 0 ][ I -B^T ][ I 0 ]
82 // [ 0 I ][ 0 I ][ 0 -\gamma S^-1 ]
83 //
84 // where the original Jacobian has the form
85 //
86 // J = [ K B^T ]
87 // [ B 0 ]
88 //
89 // and K^-1 is an approximation of the inverse of the displacement part of the
90 // Jacobian and S^-1 is an approximation of the inverse of the Schur
91 // complement S = B K^-1 B^T. The Schur complement is approximated using
92 // a mass matrix of the pressure variables.
93 class JacobianPreconditioner : public Solver
94 {
95 protected:
96  // Finite element spaces for setting up preconditioner blocks
98 
99  // Offsets for extracting block vector segments
100  Array<int> &block_trueOffsets;
101 
102  // Jacobian for block access
103  BlockOperator *jacobian;
104 
105  // Scaling factor for the pressure mass matrix in the block preconditioner
106  double gamma;
107 
108  // Objects for the block preconditioner application
109  SparseMatrix *pressure_mass;
110  Solver *mass_pcg;
111  Solver *mass_prec;
112  Solver *stiff_pcg;
113  Solver *stiff_prec;
114 
115 public:
116  JacobianPreconditioner(Array<FiniteElementSpace *> &fes,
117  SparseMatrix &mass, Array<int> &offsets);
118 
119  virtual void Mult(const Vector &k, Vector &y) const;
120  virtual void SetOperator(const Operator &op);
121 
122  virtual ~JacobianPreconditioner();
123 };
124 
125 // After spatial discretization, the rubber model can be written as:
126 // 0 = H(x)
127 // where x is the block vector representing the deformation and pressure and
128 // H(x) is the nonlinear incompressible neo-Hookean operator.
129 class RubberOperator : public Operator
130 {
131 protected:
132  // Finite element spaces
134 
135  // Block nonlinear form
136  BlockNonlinearForm *Hform;
137 
138  // Pressure mass matrix for the preconditioner
139  SparseMatrix *pressure_mass;
140 
141  // Newton solver for the hyperelastic operator
142  NewtonSolver newton_solver;
143  GeneralResidualMonitor newton_monitor;
144 
145  // Solver for the Jacobian solve in the Newton method
146  Solver *j_solver;
147  GeneralResidualMonitor j_monitor;
148 
149  // Preconditioner for the Jacobian
150  Solver *j_prec;
151 
152  // Shear modulus coefficient
153  Coefficient &mu;
154 
155  // Block offsets for variable access
156  Array<int> &block_trueOffsets;
157 
158 public:
159  RubberOperator(Array<FiniteElementSpace *> &fes, Array<Array<int> *>&ess_bdr,
160  Array<int> &block_trueOffsets, double rel_tol, double abs_tol,
161  int iter, Coefficient &mu);
162 
163  // Required to use the native newton solver
164  virtual Operator &GetGradient(const Vector &xp) const;
165  virtual void Mult(const Vector &k, Vector &y) const;
166 
167  // Driver for the newton solver
168  void Solve(Vector &xp) const;
169 
170  virtual ~RubberOperator();
171 };
172 
173 // Visualization driver
174 void visualize(ostream &out, Mesh *mesh, GridFunction *deformed_nodes,
175  GridFunction *field, const char *field_name = NULL,
176  bool init_vis = false);
177 
178 // Configuration definition functions
179 void ReferenceConfiguration(const Vector &x, Vector &y);
180 void InitialDeformation(const Vector &x, Vector &y);
181 
182 
183 int main(int argc, char *argv[])
184 {
185  // 1. Parse command-line options
186  const char *mesh_file = "../data/beam-tet.mesh";
187  int ref_levels = 0;
188  int order = 2;
189  bool visualization = true;
190  double newton_rel_tol = 1e-4;
191  double newton_abs_tol = 1e-6;
192  int newton_iter = 500;
193  double mu = 1.0;
194 
195  OptionsParser args(argc, argv);
196  args.AddOption(&mesh_file, "-m", "--mesh",
197  "Mesh file to use.");
198  args.AddOption(&ref_levels, "-r", "--refine",
199  "Number of times to refine the mesh uniformly.");
200  args.AddOption(&order, "-o", "--order",
201  "Order (degree) of the finite elements.");
202  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
203  "--no-visualization",
204  "Enable or disable GLVis visualization.");
205  args.AddOption(&newton_rel_tol, "-rel", "--relative-tolerance",
206  "Relative tolerance for the Newton solve.");
207  args.AddOption(&newton_abs_tol, "-abs", "--absolute-tolerance",
208  "Absolute tolerance for the Newton solve.");
209  args.AddOption(&newton_iter, "-it", "--newton-iterations",
210  "Maximum iterations for the Newton solve.");
211  args.AddOption(&mu, "-mu", "--shear-modulus",
212  "Shear modulus for the neo-Hookean material.");
213  args.Parse();
214  if (!args.Good())
215  {
216  args.PrintUsage(cout);
217  return 1;
218  }
219  args.PrintOptions(cout);
220 
221  // 2. Read the mesh from the given mesh file. We can handle triangular,
222  // quadrilateral, tetrahedral and hexahedral meshes with the same code.
223  Mesh *mesh = new Mesh(mesh_file, 1, 1);
224  int dim = mesh->Dimension();
225 
226  // 3. Refine the mesh to increase the resolution. In this example we do
227  // 'ref_levels' of uniform refinement, where 'ref_levels' is a
228  // command-line parameter.
229  for (int lev = 0; lev < ref_levels; lev++)
230  {
231  mesh->UniformRefinement();
232  }
233 
234  // 4. Define the shear modulus for the incompressible Neo-Hookean material
235  ConstantCoefficient c_mu(mu);
236 
237  // 5. Define the finite element spaces for displacement and pressure
238  // (Taylor-Hood elements). By default, the displacement (u/x) is a second
239  // order vector field, while the pressure (p) is a linear scalar function.
240  H1_FECollection quad_coll(order, dim);
241  H1_FECollection lin_coll(order-1, dim);
242 
243  FiniteElementSpace R_space(mesh, &quad_coll, dim, Ordering::byVDIM);
244  FiniteElementSpace W_space(mesh, &lin_coll);
245 
246  Array<FiniteElementSpace *> spaces(2);
247  spaces[0] = &R_space;
248  spaces[1] = &W_space;
249 
250  int R_size = R_space.GetTrueVSize();
251  int W_size = W_space.GetTrueVSize();
252 
253  // 6. Define the Dirichlet conditions (set to boundary attribute 1 and 2)
254  Array<Array<int> *> ess_bdr(2);
255 
256  Array<int> ess_bdr_u(R_space.GetMesh()->bdr_attributes.Max());
257  Array<int> ess_bdr_p(W_space.GetMesh()->bdr_attributes.Max());
258 
259  ess_bdr_p = 0;
260  ess_bdr_u = 0;
261  ess_bdr_u[0] = 1;
262  ess_bdr_u[1] = 1;
263 
264  ess_bdr[0] = &ess_bdr_u;
265  ess_bdr[1] = &ess_bdr_p;
266 
267  // 7. Print the mesh statistics
268  std::cout << "***********************************************************\n";
269  std::cout << "dim(u) = " << R_size << "\n";
270  std::cout << "dim(p) = " << W_size << "\n";
271  std::cout << "dim(u+p) = " << R_size + W_size << "\n";
272  std::cout << "***********************************************************\n";
273 
274  // 8. Define the block structure of the solution vector (u then p)
275  Array<int> block_trueOffsets(3);
276  block_trueOffsets[0] = 0;
277  block_trueOffsets[1] = R_space.GetTrueVSize();
278  block_trueOffsets[2] = W_space.GetTrueVSize();
279  block_trueOffsets.PartialSum();
280 
281  BlockVector xp(block_trueOffsets);
282 
283  // 9. Define grid functions for the current configuration, reference
284  // configuration, final deformation, and pressure
285  GridFunction x_gf(&R_space);
286  GridFunction x_ref(&R_space);
287  GridFunction x_def(&R_space);
288  GridFunction p_gf(&W_space);
289 
290  x_gf.MakeTRef(&R_space, xp.GetBlock(0), 0);
291  p_gf.MakeTRef(&W_space, xp.GetBlock(1), 0);
292 
295 
296  x_gf.ProjectCoefficient(deform);
297  x_ref.ProjectCoefficient(refconfig);
298  p_gf = 0.0;
299 
300  x_gf.SetTrueVector();
301  p_gf.SetTrueVector();
302 
303  // 10. Initialize the incompressible neo-Hookean operator
304  RubberOperator oper(spaces, ess_bdr, block_trueOffsets,
305  newton_rel_tol, newton_abs_tol, newton_iter, c_mu);
306 
307  // 11. Solve the Newton system
308  oper.Solve(xp);
309 
310  // 12. Compute the final deformation
311  x_gf.SetFromTrueVector();
312  p_gf.SetFromTrueVector();
313  subtract(x_gf, x_ref, x_def);
314 
315  // 13. Visualize the results if requested
316  socketstream vis_u, vis_p;
317  if (visualization)
318  {
319  char vishost[] = "localhost";
320  int visport = 19916;
321  vis_u.open(vishost, visport);
322  vis_u.precision(8);
323  visualize(vis_u, mesh, &x_gf, &x_def, "Deformation", true);
324  vis_p.open(vishost, visport);
325  vis_p.precision(8);
326  visualize(vis_p, mesh, &x_gf, &p_gf, "Pressure", true);
327  }
328 
329  // 14. Save the displaced mesh, the final deformation, and the pressure
330  {
331  GridFunction *nodes = &x_gf;
332  int owns_nodes = 0;
333  mesh->SwapNodes(nodes, owns_nodes);
334 
335  ofstream mesh_ofs("deformed.mesh");
336  mesh_ofs.precision(8);
337  mesh->Print(mesh_ofs);
338 
339  ofstream pressure_ofs("pressure.sol");
340  pressure_ofs.precision(8);
341  p_gf.Save(pressure_ofs);
342 
343  ofstream deformation_ofs("deformation.sol");
344  deformation_ofs.precision(8);
345  x_def.Save(deformation_ofs);
346  }
347 
348  // 15. Free the used memory
349  delete mesh;
350 
351  return 0;
352 }
353 
354 
355 JacobianPreconditioner::JacobianPreconditioner(Array<FiniteElementSpace *> &fes,
356  SparseMatrix &mass,
357  Array<int> &offsets)
358  : Solver(offsets[2]), block_trueOffsets(offsets), pressure_mass(&mass)
359 {
360  fes.Copy(spaces);
361 
362  gamma = 0.00001;
363 
364  // The mass matrix and preconditioner do not change every Newton cycle, so we
365  // only need to define them once
366  GSSmoother *mass_prec_gs = new GSSmoother(*pressure_mass);
367 
368  mass_prec = mass_prec_gs;
369 
370  CGSolver *mass_pcg_iter = new CGSolver();
371  mass_pcg_iter->SetRelTol(1e-12);
372  mass_pcg_iter->SetAbsTol(1e-12);
373  mass_pcg_iter->SetMaxIter(200);
374  mass_pcg_iter->SetPrintLevel(0);
375  mass_pcg_iter->SetPreconditioner(*mass_prec);
376  mass_pcg_iter->SetOperator(*pressure_mass);
377  mass_pcg_iter->iterative_mode = false;
378 
379  mass_pcg = mass_pcg_iter;
380 
381  // The stiffness matrix does change every Newton cycle, so we will define it
382  // during SetOperator
383  stiff_pcg = NULL;
384  stiff_prec = NULL;
385 }
386 
387 void JacobianPreconditioner::Mult(const Vector &k, Vector &y) const
388 {
389  // Extract the blocks from the input and output vectors
390  Vector disp_in(k.GetData() + block_trueOffsets[0],
391  block_trueOffsets[1]-block_trueOffsets[0]);
392  Vector pres_in(k.GetData() + block_trueOffsets[1],
393  block_trueOffsets[2]-block_trueOffsets[1]);
394 
395  Vector disp_out(y.GetData() + block_trueOffsets[0],
396  block_trueOffsets[1]-block_trueOffsets[0]);
397  Vector pres_out(y.GetData() + block_trueOffsets[1],
398  block_trueOffsets[2]-block_trueOffsets[1]);
399 
400  Vector temp(block_trueOffsets[1]-block_trueOffsets[0]);
401  Vector temp2(block_trueOffsets[1]-block_trueOffsets[0]);
402 
403  // Perform the block elimination for the preconditioner
404  mass_pcg->Mult(pres_in, pres_out);
405  pres_out *= -gamma;
406 
407  jacobian->GetBlock(0,1).Mult(pres_out, temp);
408  subtract(disp_in, temp, temp2);
409 
410  stiff_pcg->Mult(temp2, disp_out);
411 }
412 
413 void JacobianPreconditioner::SetOperator(const Operator &op)
414 {
415  jacobian = (BlockOperator *) &op;
416 
417  // Initialize the stiffness preconditioner and solver
418  if (stiff_prec == NULL)
419  {
420  GSSmoother *stiff_prec_gs = new GSSmoother();
421 
422  stiff_prec = stiff_prec_gs;
423 
424  GMRESSolver *stiff_pcg_iter = new GMRESSolver();
425  stiff_pcg_iter->SetRelTol(1e-8);
426  stiff_pcg_iter->SetAbsTol(1e-8);
427  stiff_pcg_iter->SetMaxIter(200);
428  stiff_pcg_iter->SetPrintLevel(0);
429  stiff_pcg_iter->SetPreconditioner(*stiff_prec);
430  stiff_pcg_iter->iterative_mode = false;
431 
432  stiff_pcg = stiff_pcg_iter;
433  }
434 
435  // At each Newton cycle, compute the new stiffness preconditioner by updating
436  // the iterative solver which, in turn, updates its preconditioner
437  stiff_pcg->SetOperator(jacobian->GetBlock(0,0));
438 }
439 
440 JacobianPreconditioner::~JacobianPreconditioner()
441 {
442  delete mass_pcg;
443  delete mass_prec;
444  delete stiff_prec;
445  delete stiff_pcg;
446 }
447 
448 
449 RubberOperator::RubberOperator(Array<FiniteElementSpace *> &fes,
450  Array<Array<int> *> &ess_bdr,
451  Array<int> &offsets,
452  double rel_tol,
453  double abs_tol,
454  int iter,
455  Coefficient &c_mu)
456  : Operator(fes[0]->GetTrueVSize() + fes[1]->GetTrueVSize()),
457  newton_solver(), newton_monitor("Newton", 1),
458  j_monitor(" GMRES", 3), mu(c_mu), block_trueOffsets(offsets)
459 {
460  Array<Vector *> rhs(2);
461  rhs = NULL; // Set all entries in the array
462 
463  fes.Copy(spaces);
464 
465  // Define the block nonlinear form
466  Hform = new BlockNonlinearForm(spaces);
467 
468  // Add the incompressible neo-Hookean integrator
469  Hform->AddDomainIntegrator(new IncompressibleNeoHookeanIntegrator(mu));
470 
471  // Set the essential boundary conditions
472  Hform->SetEssentialBC(ess_bdr, rhs);
473 
474  // Compute the pressure mass stiffness matrix
475  BilinearForm *a = new BilinearForm(spaces[1]);
476  ConstantCoefficient one(1.0);
477  a->AddDomainIntegrator(new MassIntegrator(one));
478  a->Assemble();
479  a->Finalize();
480 
481  OperatorPtr op;
482  Array<int> p_ess_tdofs;
483  a->FormSystemMatrix(p_ess_tdofs, op);
484  pressure_mass = a->LoseMat();
485  delete a;
486 
487  // Initialize the Jacobian preconditioner
488  JacobianPreconditioner *jac_prec =
489  new JacobianPreconditioner(fes, *pressure_mass, block_trueOffsets);
490  j_prec = jac_prec;
491 
492  // Set up the Jacobian solver
493  GMRESSolver *j_gmres = new GMRESSolver();
494  j_gmres->iterative_mode = false;
495  j_gmres->SetRelTol(1e-12);
496  j_gmres->SetAbsTol(1e-12);
497  j_gmres->SetMaxIter(300);
498  j_gmres->SetPrintLevel(-1);
499  j_gmres->SetMonitor(j_monitor);
500  j_gmres->SetPreconditioner(*j_prec);
501  j_solver = j_gmres;
502 
503  // Set the newton solve parameters
504  newton_solver.iterative_mode = true;
505  newton_solver.SetSolver(*j_solver);
506  newton_solver.SetOperator(*this);
507  newton_solver.SetPrintLevel(-1);
508  newton_solver.SetMonitor(newton_monitor);
509  newton_solver.SetRelTol(rel_tol);
510  newton_solver.SetAbsTol(abs_tol);
511  newton_solver.SetMaxIter(iter);
512 }
513 
514 // Solve the Newton system
515 void RubberOperator::Solve(Vector &xp) const
516 {
517  Vector zero;
518  newton_solver.Mult(zero, xp);
519  MFEM_VERIFY(newton_solver.GetConverged(),
520  "Newton Solver did not converge.");
521 }
522 
523 // compute: y = H(x,p)
524 void RubberOperator::Mult(const Vector &k, Vector &y) const
525 {
526  Hform->Mult(k, y);
527 }
528 
529 // Compute the Jacobian from the nonlinear form
530 Operator &RubberOperator::GetGradient(const Vector &xp) const
531 {
532  return Hform->GetGradient(xp);
533 }
534 
535 RubberOperator::~RubberOperator()
536 {
537  delete Hform;
538  delete pressure_mass;
539  delete j_solver;
540  delete j_prec;
541 }
542 
543 
544 // Inline visualization
545 void visualize(ostream &out, Mesh *mesh, GridFunction *deformed_nodes,
546  GridFunction *field, const char *field_name, bool init_vis)
547 {
548  if (!out)
549  {
550  return;
551  }
552 
553  GridFunction *nodes = deformed_nodes;
554  int owns_nodes = 0;
555 
556  mesh->SwapNodes(nodes, owns_nodes);
557 
558  out << "solution\n" << *mesh << *field;
559 
560  mesh->SwapNodes(nodes, owns_nodes);
561 
562  if (init_vis)
563  {
564  out << "window_size 800 800\n";
565  out << "window_title '" << field_name << "'\n";
566  if (mesh->SpaceDimension() == 2)
567  {
568  out << "view 0 0\n"; // view from top
569  out << "keys jlA\n"; // turn off perspective and light, +anti-aliasing
570  }
571  out << "keys cmA\n"; // show colorbar and mesh, +anti-aliasing
572  out << "autoscale value\n"; // update value-range; keep mesh-extents fixed
573  }
574  out << flush;
575 }
576 
578 {
579  // Set the reference, stress free, configuration
580  y = x;
581 }
582 
583 void InitialDeformation(const Vector &x, Vector &y)
584 {
585  // Set the initial configuration. Having this different from the reference
586  // configuration can help convergence
587  y = x;
588  y[1] = x[1] + 0.25*x[0];
589 }
void visualize(ostream &out, Mesh *mesh, GridFunction *deformed_nodes, GridFunction *field, const char *field_name=NULL, bool init_vis=false)
Definition: ex10.cpp:379
void InitialDeformation(const Vector &x, Vector &y)
Definition: ex10.cpp:591
virtual void Print(std::ostream &out=mfem::out) const
Definition: mesh.hpp:1387
Conjugate gradient method.
Definition: solvers.hpp:316
void ReferenceConfiguration(const Vector &x, Vector &y)
Definition: ex19.cpp:577
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
A class to handle Vectors in a block fashion.
Definition: blockvector.hpp:30
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
void SetFromTrueVector()
Shortcut for calling SetFromTrueDofs() with GetTrueVector() as argument.
Definition: gridfunc.hpp:143
A coefficient that is constant across space and time.
Definition: coefficient.hpp:78
virtual Operator & GetGradient(const Vector &x) const
Evaluate the gradient operator at the point x. The default behavior in class Operator is to generate ...
Definition: operator.hpp:98
void SwapNodes(GridFunction *&nodes, int &own_nodes_)
Definition: mesh.cpp:7386
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:476
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
void Copy(Array &copy) const
Create a copy of the internal array to the provided copy.
Definition: array.hpp:851
A class representing a general block nonlinear operator defined on the Cartesian product of multiple ...
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:652
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:199
Data type for Gauss-Seidel smoother of sparse matrix.
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:3484
void SetTrueVector()
Shortcut for calling GetTrueDofs() with GetTrueVector() as argument.
Definition: gridfunc.hpp:137
void SetPrintLevel(int print_lvl)
Definition: solvers.cpp:71
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:98
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:150
Data type sparse matrix.
Definition: sparsemat.hpp:41
constexpr char vishost[]
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:410
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9143
constexpr int visport
void MakeTRef(FiniteElementSpace *f, double *tv)
Associate a new FiniteElementSpace and new true-dof data with the GridFunction.
Definition: gridfunc.cpp:219
void SetMaxIter(int max_it)
Definition: solvers.hpp:100
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
Newton&#39;s method for solving F(x)=b for a given operator F.
Definition: solvers.hpp:464
virtual void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)
Form the linear system matrix A, see FormLinearSystem() for details.
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:543
int Dimension() const
Definition: mesh.hpp:911
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:457
A general vector function coefficient.
int SpaceDimension() const
Definition: mesh.hpp:912
Abstract base class for an iterative solver monitor.
Definition: solvers.hpp:36
void SetAbsTol(double atol)
Definition: solvers.hpp:99
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:204
void SetRelTol(double rtol)
Definition: solvers.hpp:98
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:87
SparseMatrix * LoseMat()
Nullifies the internal matrix and returns a pointer to it. Used for transfering ownership.
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:39
void subtract(const Vector &x, const Vector &y, Vector &z)
Definition: vector.cpp:438
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
Add a boolean option and set &#39;var&#39; to receive the value. Enable/disable tags are used to set the bool...
Definition: optparser.hpp:82
void PartialSum()
Fill the entries of the array with the cumulative sum of the entries.
Definition: array.cpp:103
GMRES method.
Definition: solvers.hpp:348
void SetMonitor(IterativeSolverMonitor &m)
Set the iterative solver monitor.
Definition: solvers.hpp:114
double a
Definition: lissajous.cpp:41
A &quot;square matrix&quot; operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
double mu
Definition: ex25.cpp:139
int dim
Definition: ex24.cpp:53
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:327
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: gridfunc.cpp:2278
int open(const char hostname[], int port)
Open the socket stream on &#39;port&#39; at &#39;hostname&#39;.
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:330
Vector data type.
Definition: vector.hpp:60
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:92
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:216
Base class for solvers.
Definition: operator.hpp:648
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:66
Abstract operator.
Definition: operator.hpp:24
A class to handle Block systems in a matrix-free implementation.
int main()
Vector & GetBlock(int i)
Get the i-th vector in the block.
Definition: blockvector.hpp:90
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150