MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
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
37using namespace std;
38using namespace mfem;
39
43
44class LogarithmGridFunctionCoefficient : public Coefficient
45{
46protected:
47 GridFunction *u; // grid function
48 Coefficient *obstacle;
49 real_t min_val;
50
51public:
52 LogarithmGridFunctionCoefficient(GridFunction &u_, Coefficient &obst_,
53 real_t min_val_=-36)
54 : u(&u_), obstacle(&obst_), min_val(min_val_) { }
55
56 virtual real_t Eval(ElementTransformation &T, const IntegrationPoint &ip);
57};
58
59class ExponentialGridFunctionCoefficient : public Coefficient
60{
61protected:
62 GridFunction *u;
63 Coefficient *obstacle;
64 real_t min_val;
65 real_t max_val;
66
67public:
68 ExponentialGridFunctionCoefficient(GridFunction &u_, Coefficient &obst_,
69 real_t min_val_=0.0, real_t max_val_=1e6)
70 : u(&u_), obstacle(&obst_), min_val(min_val_), max_val(max_val_) { }
71
72 virtual real_t Eval(ElementTransformation &T, const IntegrationPoint &ip);
73};
74
75int 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 real_t alpha = 1.0;
82 real_t 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 real_t 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 real_t r0 = 1.0;
163 real_t 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);
194 u_gf.ProjectCoefficient(IC_coef);
195 u_old_gf = u_gf;
196
197 // 9. Initialize the slack variable ψₕ = ln(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 real_t 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 real_t Newton_update_size = u_tmp.ComputeL2Error(zero);
304 u_tmp = u_gf;
305
306 real_t 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 real_t 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 real_t 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
383real_t LogarithmGridFunctionCoefficient::Eval(ElementTransformation &T,
384 const IntegrationPoint &ip)
385{
386 MFEM_ASSERT(u != NULL, "grid function is not set");
387
388 real_t val = u->GetValue(T, ip) - obstacle->Eval(T, ip);
389 return max(min_val, log(val));
390}
391
392real_t ExponentialGridFunctionCoefficient::Eval(ElementTransformation &T,
393 const IntegrationPoint &ip)
394{
395 MFEM_ASSERT(u != NULL, "grid function is not set");
396
397 real_t val = u->GetValue(T, ip);
398 return min(max_val, max(min_val, exp(val) + obstacle->Eval(T, ip)));
399}
400
402{
403 real_t x = pt(0), y = pt(1);
404 real_t r = sqrt(x*x + y*y);
405 real_t r0 = 0.5;
406 real_t beta = 0.9;
407
408 real_t b = r0*beta;
409 real_t tmp = sqrt(r0*r0 - b*b);
410 real_t B = tmp + b*b/tmp;
411 real_t 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 real_t x = pt(0), y = pt(1);
426 real_t r = sqrt(x*x + y*y);
427 real_t r0 = 0.5;
428 real_t a = 0.348982574111686;
429 real_t 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 real_t x = pt(0), y = pt(1);
444 real_t r = sqrt(x*x + y*y);
445 real_t r0 = 0.5;
446 real_t a = 0.348982574111686;
447 real_t 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}
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition array.hpp:697
int Size() const
Return the logical size of the array.
Definition array.hpp:144
void PartialSum()
Fill the entries of the array with the cumulative sum of the entries.
Definition array.cpp:103
A "square matrix" 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 if the AssemblyLevel is AssemblyLevel::LEGACY....
void SetDiagonalPolicy(DiagonalPolicy policy)
Sets Operator::DiagonalPolicy used upon construction of the linear system. Policies include:
Adds new Domain Integrator. Assumes ownership of bfi.
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
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.
const SparseMatrix & SpMat() const
Returns a const reference to the sparse matrix: .
A class to handle Block diagonal preconditioners in a matrix-free implementation.
void SetDiagonalBlock(int iblock, Operator *op)
Add a square block op in the block-entry (iblock, iblock).
A class to handle Block systems in a matrix-free implementation.
void SetBlock(int iRow, int iCol, Operator *op, real_t c=1.0)
Add a block op in the block-entry (iblock, jblock).
A class to handle Vectors in a block fashion.
Vector & GetBlock(int i)
Get the i-th vector in the block.
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
virtual real_t Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the coefficient in the element described by T at the point ip.
A coefficient that is constant across space and time.
Class for domain integration .
Definition lininteg.hpp:109
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:220
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition fespace.hpp:716
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition fespace.hpp:713
A general function coefficient.
Data type for Gauss-Seidel smoother of sparse matrix.
Vector coefficient defined as the Gradient of a scalar GridFunction.
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
Class for grid function - Vector with associated FE space.
Definition gridfunc.hpp:31
virtual real_t GetValue(int i, const IntegrationPoint &ip, int vdim=1) const
Definition gridfunc.cpp:449
virtual real_t ComputeH1Error(Coefficient *exsol, VectorCoefficient *exgrad, Coefficient *ell_coef, real_t Nu, int norm_type) const
virtual void MakeRef(FiniteElementSpace *f, real_t *v)
Make the GridFunction reference external data on a new FiniteElementSpace.
Definition gridfunc.cpp:203
virtual real_t ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
Class for integration point with weight.
Definition intrules.hpp:35
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition fe_coll.hpp:330
Vector with associated FE space and LinearFormIntegrators.
Adds new Domain Integrator. Assumes ownership of lfi.
void Update()
Update the object according to the associated FE space fes.
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Mesh data type.
Definition mesh.hpp:56
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition mesh.hpp:282
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
void GetNodes(Vector &node_coord) const
Definition mesh.cpp:8973
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:6211
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization if the AssemblyLevel is AssemblyLevel::LEGACY.
void Assemble(int skip_zeros=1)
void EliminateTrialDofs(const Array< int > &bdr_attr_is_ess, const Vector &sol, Vector &rhs)
Eliminate essential boundary DOFs from the columns of the system.
const SparseMatrix & SpMat() const
Returns a const reference to the sparse matrix: .
Adds a domain integrator. Assumes ownership of bfi.
@ DIAG_ONE
Set the diagonal value to one.
Definition operator.hpp:50
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void PrintUsage(std::ostream &out) const
Print the usage message.
void PrintOptions(std::ostream &out) const
Print the options.
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 'var' to receive the value. Enable/disable tags are used to set the bool...
Definition optparser.hpp:82
bool Good() const
Return true if the command line options were parsed successfully.
Scalar coefficient defined as the product of two scalar coefficients or a scalar and a scalar coeffic...
Data type sparse matrix.
Definition sparsemat.hpp:51
Scalar coefficient defined as the linear combination of two scalar coefficients or a scalar and a sca...
A general vector function coefficient.
Vector data type.
Definition vector.hpp:80
int open(const char hostname[], int port)
Open the socket stream on 'port' at 'hostname'.
Vector beta
const real_t alpha
Definition ex15.cpp:369
int dim
Definition ex24.cpp:53
Definition ex36.cpp:441
real_t spherical_obstacle(const Vector &pt)
Definition ex36.cpp:401
real_t exact_solution_obstacle(const Vector &pt)
Definition ex36.cpp:423
int main()
real_t b
Definition lissajous.cpp:42
real_t a
Definition lissajous.cpp:41
const int visport
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
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 Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
Definition table.cpp:414
int GMRES(const Operator &A, Vector &x, const Vector &b, Solver &M, int &max_iter, int m, real_t &tol, real_t atol, int printit)
GMRES method. (tolerances are squared)
Definition solvers.cpp:1342
float real_t
Definition config.hpp:43
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30
const char vishost[]