MFEM  v4.6.0 Finite element discretization library
ex36p.cpp
Go to the documentation of this file.
1 // MFEM Example 36 - Parallel Version
2 //
3 // Compile with: make ex36p
4 //
5 // Sample runs: mpirun -np 4 ex36p -o 2
6 // mpirun -np 4 ex36p -o 2 -r 4
7 //
8 // Description: This example code demonstrates the use of MFEM to solve the
9 // bound-constrained energy minimization problem
10 //
11 // minimize ||∇u||² subject to u ≥ ϕ in H¹₀.
12 //
13 // This is known as the obstacle problem, and it is a simple
14 // mathematical model for contact mechanics.
15 //
16 // In this example, the obstacle ϕ is a half-sphere centered
17 // at the origin of a circular domain Ω. After solving to a
18 // specified tolerance, the numerical solution is compared to
19 // a closed-form exact solution to assess accuracy.
20 //
21 // The problem is discretized and solved using the proximal
22 // Galerkin finite element method, introduced by Keith and
23 // Surowiec [1].
24 //
25 // This example highlights the ability of MFEM to deliver high-
26 // order solutions to variation inequality problems and
27 // showcases how to set up and solve nonlinear mixed methods.
28 //
29 // [1] Keith, B. and Surowiec, T. (2023) Proximal Galerkin: A structure-
30 // preserving finite element method for pointwise bound constraints.
31 // arXiv:2307.12444 [math.NA]
32
33 #include "mfem.hpp"
34 #include <fstream>
35 #include <iostream>
36
37 using namespace std;
38 using namespace mfem;
39
40 double spherical_obstacle(const Vector &pt);
41 double exact_solution_obstacle(const Vector &pt);
42 void exact_solution_gradient_obstacle(const Vector &pt, Vector &grad);
43
44 class LogarithmGridFunctionCoefficient : public Coefficient
45 {
46 protected:
47  GridFunction *u; // grid function
48  Coefficient *obstacle;
49  double min_val;
50
51 public:
52  LogarithmGridFunctionCoefficient(GridFunction &u_, Coefficient &obst_,
53  double min_val_=-36)
54  : u(&u_), obstacle(&obst_), min_val(min_val_) { }
55
56  virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
57 };
58
59 class ExponentialGridFunctionCoefficient : public Coefficient
60 {
61 protected:
62  GridFunction *u;
63  Coefficient *obstacle;
64  double min_val;
65  double max_val;
66
67 public:
68  ExponentialGridFunctionCoefficient(GridFunction &u_, Coefficient &obst_,
69  double min_val_=0.0, double max_val_=1e6)
70  : u(&u_), obstacle(&obst_), min_val(min_val_), max_val(max_val_) { }
71
72  virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
73 };
74
75 int main(int argc, char *argv[])
76 {
77  // 0. Initialize MPI and HYPRE.
78  Mpi::Init();
79  int num_procs = Mpi::WorldSize();
80  int myid = Mpi::WorldRank();
81  Hypre::Init();
82
83  // 1. Parse command-line options.
84  int order = 1;
85  int max_it = 10;
86  int ref_levels = 3;
87  double alpha = 1.0;
88  double tol = 1e-5;
89  bool visualization = true;
90
91  OptionsParser args(argc, argv);
92  args.AddOption(&order, "-o", "--order",
93  "Finite element order (polynomial degree).");
94  args.AddOption(&ref_levels, "-r", "--refs",
95  "Number of h-refinements.");
96  args.AddOption(&max_it, "-mi", "--max-it",
97  "Maximum number of iterations");
98  args.AddOption(&tol, "-tol", "--tol",
99  "Stopping criteria based on the difference between"
100  "successive solution updates");
101  args.AddOption(&alpha, "-step", "--step",
102  "Step size alpha");
103  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
104  "--no-visualization",
105  "Enable or disable GLVis visualization.");
106  args.Parse();
107  if (!args.Good())
108  {
109  if (myid == 0)
110  {
111  args.PrintUsage(cout);
112  }
113  return 1;
114  }
115  if (myid == 0)
116  {
117  args.PrintOptions(cout);
118  }
119
120  // 2. Read the mesh from the mesh file.
121  const char *mesh_file = "../data/disc-nurbs.mesh";
122  Mesh mesh(mesh_file, 1, 1);
123  int dim = mesh.Dimension();
124
125  // 3. Postprocess the mesh.
126  // 3A. Refine the mesh to increase the resolution.
127  for (int l = 0; l < ref_levels; l++)
128  {
129  mesh.UniformRefinement();
130  }
131
132  // 3B. Interpolate the geometry after refinement to control geometry error.
133  // NOTE: Minimum second-order interpolation is used to improve the accuracy.
134  int curvature_order = max(order,2);
135  mesh.SetCurvature(curvature_order);
136
137  // 3C. Rescale the domain to a unit circle (radius = 1).
138  GridFunction *nodes = mesh.GetNodes();
139  double scale = 2*sqrt(2);
140  *nodes /= scale;
141
142  ParMesh pmesh(MPI_COMM_WORLD, mesh);
143  mesh.Clear();
144
145  // 4. Define the necessary finite element spaces on the mesh.
146  H1_FECollection H1fec(order+1, dim);
147  ParFiniteElementSpace H1fes(&pmesh, &H1fec);
148
149  L2_FECollection L2fec(order-1, dim);
150  ParFiniteElementSpace L2fes(&pmesh, &L2fec);
151
152  int num_dofs_H1 = H1fes.GetTrueVSize();
153  MPI_Allreduce(MPI_IN_PLACE, &num_dofs_H1, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
154  int num_dofs_L2 = L2fes.GetTrueVSize();
155  MPI_Allreduce(MPI_IN_PLACE, &num_dofs_L2, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
156  if (myid == 0)
157  {
158  cout << "Number of H1 finite element unknowns: "
159  << num_dofs_H1 << endl;
160  cout << "Number of L2 finite element unknowns: "
161  << num_dofs_L2 << endl;
162  }
163
164  Array<int> offsets(3);
165  offsets[0] = 0;
166  offsets[1] = H1fes.GetVSize();
167  offsets[2] = L2fes.GetVSize();
168  offsets.PartialSum();
169
170  Array<int> toffsets(3);
171  toffsets[0] = 0;
172  toffsets[1] = H1fes.GetTrueVSize();
173  toffsets[2] = L2fes.GetTrueVSize();
174  toffsets.PartialSum();
175
176  BlockVector x(offsets), rhs(offsets);
177  x = 0.0; rhs = 0.0;
178
179  BlockVector tx(toffsets), trhs(toffsets);
180  tx = 0.0; trhs = 0.0;
181
182  // 5. Determine the list of true (i.e. conforming) essential boundary dofs.
183  Array<int> empty;
184  Array<int> ess_tdof_list;
185  if (pmesh.bdr_attributes.Size())
186  {
187  Array<int> ess_bdr(pmesh.bdr_attributes.Max());
188  ess_bdr = 1;
189  H1fes.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
190  }
191
192  // 6. Define an initial guess for the solution.
193  auto IC_func = [](const Vector &x)
194  {
195  double r0 = 1.0;
196  double rr = 0.0;
197  for (int i=0; i<x.Size(); i++)
198  {
199  rr += x(i)*x(i);
200  }
201  return r0*r0 - rr;
202  };
203  ConstantCoefficient one(1.0);
204  ConstantCoefficient zero(0.0);
205
206  // 7. Define the solution vectors as a finite element grid functions
207  // corresponding to the fespaces.
208  ParGridFunction u_gf, delta_psi_gf;
209  u_gf.MakeRef(&H1fes,x,offsets[0]);
210  delta_psi_gf.MakeRef(&L2fes,x,offsets[1]);
211  delta_psi_gf = 0.0;
212
213  ParGridFunction u_old_gf(&H1fes);
214  ParGridFunction psi_old_gf(&L2fes);
215  ParGridFunction psi_gf(&L2fes);
216  u_old_gf = 0.0;
217  psi_old_gf = 0.0;
218
219  // 8. Define the function coefficients for the solution and use them to
220  // initialize the initial guess
223  FunctionCoefficient IC_coef(IC_func);
224  ConstantCoefficient f(0.0);
226  u_gf.ProjectCoefficient(IC_coef);
227  u_old_gf = u_gf;
228
229  // 9. Initialize the slack variable ψₕ = exp(uₕ)
230  LogarithmGridFunctionCoefficient ln_u(u_gf, obstacle);
231  psi_gf.ProjectCoefficient(ln_u);
232  psi_old_gf = psi_gf;
233
234  char vishost[] = "localhost";
235  int visport = 19916;
236  socketstream sol_sock;
237  if (visualization)
238  {
239  sol_sock.open(vishost,visport);
240  sol_sock.precision(8);
241  }
242
243  // 10. Iterate
244  int k;
245  int total_iterations = 0;
246  double increment_u = 0.1;
247  for (k = 0; k < max_it; k++)
248  {
249  ParGridFunction u_tmp(&H1fes);
250  u_tmp = u_old_gf;
251
252  if (myid == 0)
253  {
254  mfem::out << "\nOUTER ITERATION " << k+1 << endl;
255  }
256
257  int j;
258  for ( j = 0; j < 10; j++)
259  {
260  total_iterations++;
261
262  ConstantCoefficient alpha_cf(alpha);
263
264  ParLinearForm b0,b1;
265  b0.Update(&H1fes,rhs.GetBlock(0),0);
266  b1.Update(&L2fes,rhs.GetBlock(1),0);
267
268  ExponentialGridFunctionCoefficient exp_psi(psi_gf, zero);
269  ProductCoefficient neg_exp_psi(-1.0,exp_psi);
270  GradientGridFunctionCoefficient grad_u_old(&u_old_gf);
271  ProductCoefficient alpha_f(alpha, f);
272  GridFunctionCoefficient psi_cf(&psi_gf);
273  GridFunctionCoefficient psi_old_cf(&psi_old_gf);
274  SumCoefficient psi_old_minus_psi(psi_old_cf, psi_cf, 1.0, -1.0);
275
276  b0.AddDomainIntegrator(new DomainLFIntegrator(alpha_f));
277  b0.AddDomainIntegrator(new DomainLFIntegrator(psi_old_minus_psi));
278  b0.Assemble();
279
280  b1.AddDomainIntegrator(new DomainLFIntegrator(exp_psi));
281  b1.AddDomainIntegrator(new DomainLFIntegrator(obstacle));
282  b1.Assemble();
283
284  ParBilinearForm a00(&H1fes);
286  a00.AddDomainIntegrator(new DiffusionIntegrator(alpha_cf));
287  a00.Assemble();
288  HypreParMatrix A00;
289  a00.FormLinearSystem(ess_tdof_list, x.GetBlock(0), rhs.GetBlock(0),
290  A00, tx.GetBlock(0), trhs.GetBlock(0));
291
292
293  ParMixedBilinearForm a10(&H1fes,&L2fes);
295  a10.Assemble();
296  HypreParMatrix A10;
297  a10.FormRectangularLinearSystem(ess_tdof_list, empty, x.GetBlock(0),
298  rhs.GetBlock(1),
299  A10, tx.GetBlock(0), trhs.GetBlock(1));
300
301  HypreParMatrix *A01 = A10.Transpose();
302
303  ParBilinearForm a11(&L2fes);
304  a11.AddDomainIntegrator(new MassIntegrator(neg_exp_psi));
305  // NOTE: Shift the spectrum of the Hessian matrix for additional
306  // stability (Quasi-Newton).
307  ConstantCoefficient eps_cf(-1e-6);
308  if (order == 1)
309  {
310  // NOTE: ∇ₕuₕ = 0 for constant functions.
311  // Therefore, we use the mass matrix to shift the spectrum
312  a11.AddDomainIntegrator(new MassIntegrator(eps_cf));
313  }
314  else
315  {
316  a11.AddDomainIntegrator(new DiffusionIntegrator(eps_cf));
317  }
318  a11.Assemble();
319  a11.Finalize();
320  HypreParMatrix A11;
321  a11.FormSystemMatrix(empty, A11);
322
323  BlockOperator A(toffsets);
324  A.SetBlock(0,0,&A00);
325  A.SetBlock(1,0,&A10);
326  A.SetBlock(0,1,A01);
327  A.SetBlock(1,1,&A11);
328
329  BlockDiagonalPreconditioner prec(toffsets);
330  HypreBoomerAMG P00(A00);
331  P00.SetPrintLevel(0);
332  HypreSmoother P11(A11);
333  prec.SetDiagonalBlock(0,&P00);
334  prec.SetDiagonalBlock(1,&P11);
335
336  GMRESSolver gmres(MPI_COMM_WORLD);
337  gmres.SetPrintLevel(-1);
338  gmres.SetRelTol(1e-8);
339  gmres.SetMaxIter(20000);
340  gmres.SetKDim(500);
341  gmres.SetOperator(A);
342  gmres.SetPreconditioner(prec);
343  gmres.Mult(trhs,tx);
344
345  u_gf.SetFromTrueDofs(tx.GetBlock(0));
346  delta_psi_gf.SetFromTrueDofs(tx.GetBlock(1));
347
348  u_tmp -= u_gf;
349  double Newton_update_size = u_tmp.ComputeL2Error(zero);
350  u_tmp = u_gf;
351
352  double gamma = 1.0;
353  delta_psi_gf *= gamma;
354  psi_gf += delta_psi_gf;
355
356  if (visualization)
357  {
358  sol_sock << "parallel " << num_procs << " " << myid << "\n";
359  sol_sock << "solution\n" << pmesh << u_gf << "window_title 'Discrete solution'"
360  << flush;
361  }
362
363  if (myid == 0)
364  {
365  mfem::out << "Newton_update_size = " << Newton_update_size << endl;
366  }
367
368  delete A01;
369
370  if (Newton_update_size < increment_u)
371  {
372  break;
373  }
374  }
375
376  u_tmp = u_gf;
377  u_tmp -= u_old_gf;
378  increment_u = u_tmp.ComputeL2Error(zero);
379
380  if (myid == 0)
381  {
382  mfem::out << "Number of Newton iterations = " << j+1 << endl;
383  mfem::out << "Increment (|| uₕ - uₕ_prvs||) = " << increment_u << endl;
384  }
385
386  u_old_gf = u_gf;
387  psi_old_gf = psi_gf;
388
389  if (increment_u < tol || k == max_it-1)
390  {
391  break;
392  }
393
394  double H1_error = u_gf.ComputeH1Error(&exact_coef,&exact_grad_coef);
395  if (myid == 0)
396  {
397  mfem::out << "H1-error (|| u - uₕᵏ||) = " << H1_error << endl;
398  }
399
400  }
401
402  if (myid == 0)
403  {
404  mfem::out << "\n Outer iterations: " << k+1
405  << "\n Total iterations: " << total_iterations
406  << "\n Total dofs: " << num_dofs_H1 + num_dofs_L2
407  << endl;
408  }
409
410  // 11. Exact solution.
411  if (visualization)
412  {
413  socketstream err_sock(vishost, visport);
414  err_sock.precision(8);
415
416  ParGridFunction error_gf(&H1fes);
417  error_gf.ProjectCoefficient(exact_coef);
418  error_gf -= u_gf;
419
420  err_sock << "parallel " << num_procs << " " << myid << "\n";
421  err_sock << "solution\n" << pmesh << error_gf << "window_title 'Error'" <<
422  flush;
423  }
424
425  {
426  double L2_error = u_gf.ComputeL2Error(exact_coef);
427  double H1_error = u_gf.ComputeH1Error(&exact_coef,&exact_grad_coef);
428
429  ExponentialGridFunctionCoefficient u_alt_cf(psi_gf,obstacle);
430  ParGridFunction u_alt_gf(&L2fes);
431  u_alt_gf.ProjectCoefficient(u_alt_cf);
432  double L2_error_alt = u_alt_gf.ComputeL2Error(exact_coef);
433
434  if (myid == 0)
435  {
436  mfem::out << "\n Final L2-error (|| u - uₕ||) = " << L2_error <<
437  endl;
438  mfem::out << " Final H1-error (|| u - uₕ||) = " << H1_error << endl;
439  mfem::out << " Final L2-error (|| u - ϕ - exp(ψₕ)||) = " << L2_error_alt <<
440  endl;
441  }
442  }
443
444  return 0;
445 }
446
448  const IntegrationPoint &ip)
449 {
450  MFEM_ASSERT(u != NULL, "grid function is not set");
451
452  double val = u->GetValue(T, ip) - obstacle->Eval(T, ip);
453  return max(min_val, log(val));
454 }
455
457  const IntegrationPoint &ip)
458 {
459  MFEM_ASSERT(u != NULL, "grid function is not set");
460
461  double val = u->GetValue(T, ip);
462  return min(max_val, max(min_val, exp(val) + obstacle->Eval(T, ip)));
463 }
464
465 double spherical_obstacle(const Vector &pt)
466 {
467  double x = pt(0), y = pt(1);
468  double r = sqrt(x*x + y*y);
469  double r0 = 0.5;
470  double beta = 0.9;
471
472  double b = r0*beta;
473  double tmp = sqrt(r0*r0 - b*b);
474  double B = tmp + b*b/tmp;
475  double C = -b/tmp;
476
477  if (r > b)
478  {
479  return B + r * C;
480  }
481  else
482  {
483  return sqrt(r0*r0 - r*r);
484  }
485 }
486
488 {
489  double x = pt(0), y = pt(1);
490  double r = sqrt(x*x + y*y);
491  double r0 = 0.5;
492  double a = 0.348982574111686;
493  double A = -0.340129705945858;
494
495  if (r > a)
496  {
497  return A * log(r);
498  }
499  else
500  {
501  return sqrt(r0*r0-r*r);
502  }
503 }
504
506 {
507  double x = pt(0), y = pt(1);
508  double r = sqrt(x*x + y*y);
509  double r0 = 0.5;
510  double a = 0.348982574111686;
511  double A = -0.340129705945858;
512
513  if (r > a)
514  {
515  grad(0) = A * x / (r*r);
516  grad(1) = A * y / (r*r);
517  }
518  else
519  {
520  grad(0) = - x / sqrt( r0*r0 - r*r );
521  grad(1) = - y / sqrt( r0*r0 - r*r );
522  }
523 }
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:108
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Definition: pfespace.cpp:1031
double Eval(ElementTransformation &T, const IntegrationPoint &ip, double t)
Evaluate the coefficient in the element described by T at the point ip at time t. ...
Definition: coefficient.hpp:68
int visport
int main(int argc, char *argv[])
Definition: ex36p.cpp:75
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
A coefficient that is constant across space and time.
Definition: coefficient.hpp:84
virtual void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B, int copy_interior=0)
Form the linear system A X = B, corresponding to this bilinear form and the linear form b(...
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:331
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1020
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:462
Scalar coefficient defined as the product of two scalar coefficients or a scalar and a scalar coeffic...
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:980
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:68
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:159
Abstract parallel finite element space.
Definition: pfespace.hpp:28
Vector beta
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: pgridfunc.cpp:561
STL namespace.
void Update(ParFiniteElementSpace *pf=NULL)
Update the object according to the given new FE space *pf.
Definition: plinearform.cpp:21
void SetDiagonalPolicy(DiagonalPolicy policy)
Sets diagonal policy used upon construction of the linear system.
virtual void SetOperator(const Operator &op) override
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:179
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds a domain integrator. Assumes ownership of bfi.
virtual void MakeRef(FiniteElementSpace *f, double *v)
Make the ParGridFunction reference external data on a new FiniteElementSpace.
Definition: pgridfunc.cpp:111
double spherical_obstacle(const Vector &pt)
Definition: ex36p.cpp:465
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1590
Class for parallel linear form.
Definition: plinearform.hpp:26
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition: solvers.cpp:71
virtual void SetFromTrueDofs(const Vector &tv)
Set the GridFunction from the given true-dof vector.
Definition: pgridfunc.hpp:169
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
A class to handle Block diagonal preconditioners in a matrix-free implementation. ...
char vishost[]
void Assemble()
Assembles the ParLinearForm i.e. sums over all domain/bdr integrators.
Definition: plinearform.cpp:46
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:10232
void SetPrintLevel(int print_level)
Definition: hypre.hpp:1670
void SetMaxIter(int max_it)
Definition: solvers.hpp:201
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Set the curvature of the mesh nodes using the given polynomial degree.
Definition: mesh.cpp:5635
Vector coefficient defined as the Gradient of a scalar GridFunction.
void Assemble(int skip_zeros=1)
void Assemble(int skip_zeros=1)
Assemble the local matrix.
void SetKDim(int dim)
Set the number of iteration to perform between restarts, default is 50.
Definition: solvers.hpp:538
Parallel smoothers in hypre.
Definition: hypre.hpp:971
virtual double ComputeH1Error(Coefficient *exsol, VectorCoefficient *exgrad, Coefficient *ell_coef, double Nu, int norm_type) const
Definition: pgridfunc.hpp:343
Set the diagonal value to one.
Definition: operator.hpp:50
virtual double ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const
Definition: pgridfunc.hpp:285
A general vector function coefficient.
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:41
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:275
void SetRelTol(double rtol)
Definition: solvers.hpp:199
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:41
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
double exact_solution_obstacle(const Vector &pt)
Definition: ex36p.cpp:487
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:525
void exact_solution_gradient_obstacle(const Vector &pt, Vector &grad)
Definition: ex36p.cpp:505
HypreParMatrix * Transpose() const
Returns the transpose of *this.
Definition: hypre.cpp:1611
virtual int GetTrueVSize() const
Return the number of local vector true dofs.
Definition: pfespace.hpp:285
double a
Definition: lissajous.cpp:41
virtual void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)
Form the linear system matrix A, see FormLinearSystem() for details.
Class for integration point with weight.
Definition: intrules.hpp:31
Class for parallel bilinear form using different test and trial FE spaces.
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
int dim
Definition: ex24.cpp:53
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:709
virtual void FormRectangularLinearSystem(const Array< int > &trial_tdof_list, const Array< int > &test_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B)
Form the parallel linear system A X = B, corresponding to this mixed bilinear form and the linear for...
Class for parallel bilinear form.
int open(const char hostname[], int port)
Open the socket stream on &#39;port&#39; at &#39;hostname&#39;.
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
void Clear()
Clear the contents of the Mesh.
Definition: mesh.hpp:678
const double alpha
Definition: ex15.cpp:369
A general function coefficient.
Vector data type.
Definition: vector.hpp:58
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:173
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:259
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:8302
double u(const Vector &xvec)
Definition: lor_mms.hpp:22
Class for parallel grid function.
Definition: pgridfunc.hpp:32
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
A class to handle Block systems in a matrix-free implementation.
Class for parallel meshes.
Definition: pmesh.hpp:32
Scalar coefficient defined as the linear combination of two scalar coefficients or a scalar and a sca...
void SetBlock(int iRow, int iCol, Operator *op, double c=1.0)
Add a block op in the block-entry (iblock, jblock).
void SetDiagonalBlock(int iblock, Operator *op)
Add a square block op in the block-entry (iblock, iblock).
Vector & GetBlock(int i)
Get the i-th vector in the block.
Definition: blockvector.hpp:93
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:327
double f(const Vector &p)