MFEM  v4.6.0 Finite element discretization library
ex36.cpp
Go to the documentation of this file.
1 // MFEM Example 36
2 //
3 // Compile with: make ex36
4 //
5 // Sample runs: ex36 -o 2
6 // ex36 -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);
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  // 1. Parse command-line options.
78  int order = 1;
79  int max_it = 10;
80  int ref_levels = 3;
81  double alpha = 1.0;
82  double tol = 1e-5;
83  bool visualization = true;
84
85  OptionsParser args(argc, argv);
87  "Finite element order (polynomial degree).");
89  "Number of h-refinements.");
91  "Maximum number of iterations");
93  "Stopping criteria based on the difference between"
96  "Step size alpha");
98  "--no-visualization",
99  "Enable or disable GLVis visualization.");
100  args.Parse();
101  if (!args.Good())
102  {
103  args.PrintUsage(cout);
104  return 1;
105  }
106  args.PrintOptions(cout);
107
108  // 2. Read the mesh from the mesh file.
109  const char *mesh_file = "../data/disc-nurbs.mesh";
110  Mesh mesh(mesh_file, 1, 1);
111  int dim = mesh.Dimension();
112
113  // 3. Postprocess the mesh.
114  // 3A. Refine the mesh to increase the resolution.
115  for (int l = 0; l < ref_levels; l++)
116  {
117  mesh.UniformRefinement();
118  }
119
120  // 3B. Interpolate the geometry after refinement to control geometry error.
121  // NOTE: Minimum second-order interpolation is used to improve the accuracy.
122  int curvature_order = max(order,2);
123  mesh.SetCurvature(curvature_order);
124
125  // 3C. Rescale the domain to a unit circle (radius = 1).
126  GridFunction *nodes = mesh.GetNodes();
127  double scale = 2*sqrt(2);
128  *nodes /= scale;
129
130  // 4. Define the necessary finite element spaces on the mesh.
131  H1_FECollection H1fec(order+1, dim);
132  FiniteElementSpace H1fes(&mesh, &H1fec);
133
134  L2_FECollection L2fec(order-1, dim);
135  FiniteElementSpace L2fes(&mesh, &L2fec);
136
137  cout << "Number of H1 finite element unknowns: "
138  << H1fes.GetTrueVSize() << endl;
139  cout << "Number of L2 finite element unknowns: "
140  << L2fes.GetTrueVSize() << endl;
141
142  Array<int> offsets(3);
143  offsets[0] = 0;
144  offsets[1] = H1fes.GetVSize();
145  offsets[2] = L2fes.GetVSize();
146  offsets.PartialSum();
147
148  BlockVector x(offsets), rhs(offsets);
149  x = 0.0; rhs = 0.0;
150
151  // 5. Determine the list of true (i.e. conforming) essential boundary dofs.
152  Array<int> ess_bdr;
153  if (mesh.bdr_attributes.Size())
154  {
155  ess_bdr.SetSize(mesh.bdr_attributes.Max());
156  ess_bdr = 1;
157  }
158
159  // 6. Define an initial guess for the solution.
160  auto IC_func = [](const Vector &x)
161  {
162  double r0 = 1.0;
163  double rr = 0.0;
164  for (int i=0; i<x.Size(); i++)
165  {
166  rr += x(i)*x(i);
167  }
168  return r0*r0 - rr;
169  };
170  ConstantCoefficient one(1.0);
171  ConstantCoefficient zero(0.0);
172
173  // 7. Define the solution vectors as a finite element grid functions
174  // corresponding to the fespaces.
175  GridFunction u_gf, delta_psi_gf;
176
177  u_gf.MakeRef(&H1fes,x,offsets[0]);
178  delta_psi_gf.MakeRef(&L2fes,x,offsets[1]);
179  delta_psi_gf = 0.0;
180
181  GridFunction u_old_gf(&H1fes);
182  GridFunction psi_old_gf(&L2fes);
183  GridFunction psi_gf(&L2fes);
184  u_old_gf = 0.0;
185  psi_old_gf = 0.0;
186
187  // 8. Define the function coefficients for the solution and use them to
188  // initialize the initial guess
191  FunctionCoefficient IC_coef(IC_func);
192  ConstantCoefficient f(0.0);
194  u_gf.ProjectCoefficient(IC_coef);
195  u_old_gf = u_gf;
196
197  // 9. Initialize the slack variable ψₕ = exp(uₕ)
198  LogarithmGridFunctionCoefficient ln_u(u_gf, obstacle);
199  psi_gf.ProjectCoefficient(ln_u);
200  psi_old_gf = psi_gf;
201
202  char vishost[] = "localhost";
203  int visport = 19916;
204  socketstream sol_sock;
205  if (visualization)
206  {
207  sol_sock.open(vishost,visport);
208  sol_sock.precision(8);
209  }
210
211  // 10. Iterate
212  int k;
213  int total_iterations = 0;
214  double increment_u = 0.1;
215  for (k = 0; k < max_it; k++)
216  {
217  GridFunction u_tmp(&H1fes);
218  u_tmp = u_old_gf;
219
220  mfem::out << "\nOUTER ITERATION " << k+1 << endl;
221
222  int j;
223  for ( j = 0; j < 10; j++)
224  {
225  total_iterations++;
226
227  ConstantCoefficient alpha_cf(alpha);
228
229  LinearForm b0,b1;
230  b0.Update(&H1fes,rhs.GetBlock(0),0);
231  b1.Update(&L2fes,rhs.GetBlock(1),0);
232
233  ExponentialGridFunctionCoefficient exp_psi(psi_gf, zero);
234  ProductCoefficient neg_exp_psi(-1.0,exp_psi);
236  ProductCoefficient alpha_f(alpha, f);
237  GridFunctionCoefficient psi_cf(&psi_gf);
238  GridFunctionCoefficient psi_old_cf(&psi_old_gf);
239  SumCoefficient psi_old_minus_psi(psi_old_cf, psi_cf, 1.0, -1.0);
240
243  b0.Assemble();
244
247  b1.Assemble();
248
249  BilinearForm a00(&H1fes);
252  a00.Assemble();
253  a00.EliminateEssentialBC(ess_bdr,x.GetBlock(0),rhs.GetBlock(0),
255  a00.Finalize();
256  SparseMatrix &A00 = a00.SpMat();
257
258  MixedBilinearForm a10(&H1fes,&L2fes);
260  a10.Assemble();
261  a10.EliminateTrialDofs(ess_bdr, x.GetBlock(0), rhs.GetBlock(1));
262  a10.Finalize();
263  SparseMatrix &A10 = a10.SpMat();
264
265  SparseMatrix *A01 = Transpose(A10);
266
267  BilinearForm a11(&L2fes);
269  // NOTE: Shift the spectrum of the Hessian matrix for additional
270  // stability (Quasi-Newton).
271  ConstantCoefficient eps_cf(-1e-6);
272  if (order == 1)
273  {
274  // NOTE: ∇ₕuₕ = 0 for constant functions.
275  // Therefore, we use the mass matrix to shift the spectrum
277  }
278  else
279  {
281  }
282  a11.Assemble();
283  a11.Finalize();
284  SparseMatrix &A11 = a11.SpMat();
285
286  BlockOperator A(offsets);
287  A.SetBlock(0,0,&A00);
288  A.SetBlock(1,0,&A10);
289  A.SetBlock(0,1,A01);
290  A.SetBlock(1,1,&A11);
291
292  BlockDiagonalPreconditioner prec(offsets);
293  prec.SetDiagonalBlock(0,new GSSmoother(A00));
294  prec.SetDiagonalBlock(1,new GSSmoother(A11));
295  prec.owns_blocks = 1;
296
297  GMRES(A,prec,rhs,x,0,10000,500,1e-12,0.0);
298
299  u_gf.MakeRef(&H1fes, x.GetBlock(0), 0);
300  delta_psi_gf.MakeRef(&L2fes, x.GetBlock(1), 0);
301
302  u_tmp -= u_gf;
303  double Newton_update_size = u_tmp.ComputeL2Error(zero);
304  u_tmp = u_gf;
305
306  double gamma = 1.0;
307  delta_psi_gf *= gamma;
308  psi_gf += delta_psi_gf;
309
310  if (visualization)
311  {
312  sol_sock << "solution\n" << mesh << u_gf << "window_title 'Discrete solution'"
313  << flush;
314  mfem::out << "Newton_update_size = " << Newton_update_size << endl;
315  }
316
317  delete A01;
318
319  if (Newton_update_size < increment_u)
320  {
321  break;
322  }
323  }
324
325  u_tmp = u_gf;
326  u_tmp -= u_old_gf;
327  increment_u = u_tmp.ComputeL2Error(zero);
328
329  mfem::out << "Number of Newton iterations = " << j+1 << endl;
330  mfem::out << "Increment (|| uₕ - uₕ_prvs||) = " << increment_u << endl;
331
332  u_old_gf = u_gf;
333  psi_old_gf = psi_gf;
334
335  if (increment_u < tol || k == max_it-1)
336  {
337  break;
338  }
339
341  mfem::out << "H1-error (|| u - uₕᵏ||) = " << H1_error << endl;
342
343  }
344
345  mfem::out << "\n Outer iterations: " << k+1
346  << "\n Total iterations: " << total_iterations
347  << "\n Total dofs: " << H1fes.GetTrueVSize() + L2fes.GetTrueVSize()
348  << endl;
349
350  // 11. Exact solution.
351  if (visualization)
352  {
353  socketstream err_sock(vishost, visport);
354  err_sock.precision(8);
355
356  GridFunction error_gf(&H1fes);
357  error_gf.ProjectCoefficient(exact_coef);
358  error_gf -= u_gf;
359
360  err_sock << "solution\n" << mesh << error_gf << "window_title 'Error'" <<
361  flush;
362  }
363
364  {
365  double L2_error = u_gf.ComputeL2Error(exact_coef);
367
368  ExponentialGridFunctionCoefficient u_alt_cf(psi_gf,obstacle);
369  GridFunction u_alt_gf(&L2fes);
370  u_alt_gf.ProjectCoefficient(u_alt_cf);
371  double L2_error_alt = u_alt_gf.ComputeL2Error(exact_coef);
372
373  mfem::out << "\n Final L2-error (|| u - uₕ||) = " << L2_error <<
374  endl;
375  mfem::out << " Final H1-error (|| u - uₕ||) = " << H1_error << endl;
376  mfem::out << " Final L2-error (|| u - ϕ - exp(ψₕ)||) = " << L2_error_alt <<
377  endl;
378  }
379
380  return 0;
381 }
382
384  const IntegrationPoint &ip)
385 {
386  MFEM_ASSERT(u != NULL, "grid function is not set");
387
388  double val = u->GetValue(T, ip) - obstacle->Eval(T, ip);
389  return max(min_val, log(val));
390 }
391
393  const IntegrationPoint &ip)
394 {
395  MFEM_ASSERT(u != NULL, "grid function is not set");
396
397  double val = u->GetValue(T, ip);
398  return min(max_val, max(min_val, exp(val) + obstacle->Eval(T, ip)));
399 }
400
401 double spherical_obstacle(const Vector &pt)
402 {
403  double x = pt(0), y = pt(1);
404  double r = sqrt(x*x + y*y);
405  double r0 = 0.5;
406  double beta = 0.9;
407
408  double b = r0*beta;
409  double tmp = sqrt(r0*r0 - b*b);
410  double B = tmp + b*b/tmp;
411  double C = -b/tmp;
412
413  if (r > b)
414  {
415  return B + r * C;
416  }
417  else
418  {
419  return sqrt(r0*r0 - r*r);
420  }
421 }
422
424 {
425  double x = pt(0), y = pt(1);
426  double r = sqrt(x*x + y*y);
427  double r0 = 0.5;
428  double a = 0.348982574111686;
429  double A = -0.340129705945858;
430
431  if (r > a)
432  {
433  return A * log(r);
434  }
435  else
436  {
437  return sqrt(r0*r0-r*r);
438  }
439 }
440
442 {
443  double x = pt(0), y = pt(1);
444  double r = sqrt(x*x + y*y);
445  double r0 = 0.5;
446  double a = 0.348982574111686;
447  double A = -0.340129705945858;
448
449  if (r > a)
450  {
451  grad(0) = A * x / (r*r);
452  grad(1) = A * y / (r*r);
453  }
454  else
455  {
456  grad(0) = - x / sqrt( r0*r0 - r*r );
457  grad(1) = - y / sqrt( r0*r0 - r*r );
458  }
459 }
virtual double ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const
Definition: gridfunc.cpp:2786
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:108
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
double exact_solution_obstacle(const Vector &pt)
Definition: ex36.cpp:423
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.
A coefficient that is constant across space and time.
Definition: coefficient.hpp:84
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:331
const SparseMatrix & SpMat() const
Returns a const reference to the sparse matrix: .
void EliminateTrialDofs(const Array< int > &bdr_attr_is_ess, const Vector &sol, Vector &rhs)
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
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:168
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.
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
Vector beta
STL namespace.
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
double spherical_obstacle(const Vector &pt)
Definition: ex36.cpp:401
void SetDiagonalPolicy(DiagonalPolicy policy)
Sets diagonal policy used upon construction of the linear system.
Data type for Gauss-Seidel smoother of sparse matrix.
Adds a domain integrator. Assumes ownership of bfi.
const SparseMatrix & SpMat() const
Returns a const reference to the sparse matrix: .
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. ...
Data type sparse matrix.
Definition: sparsemat.hpp:50
char vishost[]
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:10232
void EliminateEssentialBC(const Array< int > &bdr_attr_is_ess, const Vector &sol, Vector &rhs, DiagonalPolicy dpolicy=DIAG_ONE)
Eliminate essential boundary DOFs from the system.
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)
int main(int argc, char *argv[])
Definition: ex36.cpp:75
virtual void MakeRef(FiniteElementSpace *f, double *v)
Make the GridFunction reference external data on a new FiniteElementSpace.
Definition: gridfunc.cpp:203
Set the diagonal value to one.
Definition: operator.hpp:50
A general vector function coefficient.
Adds new Domain Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:41
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:712
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:275
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:219
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:41
virtual double ComputeH1Error(Coefficient *exsol, VectorCoefficient *exgrad, Coefficient *ell_coef, double Nu, int norm_type) const
Definition: gridfunc.cpp:3160
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
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 SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:687
void PartialSum()
Fill the entries of the array with the cumulative sum of the entries.
Definition: array.cpp:103
void Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
Definition: table.cpp:413
void Update()
Update the object according to the associated FE space fes.
Definition: linearform.cpp:344
double a
Definition: lissajous.cpp:41
A "square matrix" operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
Class for integration point with weight.
Definition: intrules.hpp:31
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
int dim
Definition: ex24.cpp:53
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 ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: gridfunc.cpp:2415
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
const double alpha
Definition: ex15.cpp:369
A general function coefficient.
Vector data type.
Definition: vector.hpp:58
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:259
Vector with associated FE space and LinearFormIntegrators.
Definition: linearform.hpp:24
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:8302
double u(const Vector &xvec)
Definition: lor_mms.hpp:22
Definition: ex36.cpp:441
A class to handle Block systems in a matrix-free implementation.
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).
int GMRES(const Operator &A, Vector &x, const Vector &b, Solver &M, int &max_iter, int m, double &tol, double atol, int printit)
GMRES method. (tolerances are squared)
Definition: solvers.cpp:1340
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)