MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
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
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 // 0. Initialize MPI and HYPRE.
78 Mpi::Init();
79 int num_procs = Mpi::WorldSize();
80 int myid = Mpi::WorldRank();
82
83 // 1. Parse command-line options.
84 int order = 1;
85 int max_it = 10;
86 int ref_levels = 3;
87 real_t alpha = 1.0;
88 real_t tol = 1e-5;
89 bool visualization = true;
90
91 OptionsParser args(argc, argv);
93 "Finite element order (polynomial degree).");
95 "Number of h-refinements.");
97 "Maximum number of iterations");
99 "Stopping criteria based on the difference between"
102 "Step size alpha");
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 real_t 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 real_t r0 = 1.0;
196 real_t 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);
226 u_gf.ProjectCoefficient(IC_coef);
227 u_old_gf = u_gf;
228
229 // 9. Initialize the slack variable ψₕ = ln(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 real_t 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);
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
278 b0.Assemble();
279
282 b1.Assemble();
283
284 ParBilinearForm a00(&H1fes);
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);
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
313 }
314 else
315 {
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 real_t Newton_update_size = u_tmp.ComputeL2Error(zero);
350 u_tmp = u_gf;
351
352 real_t 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
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 real_t L2_error = u_gf.ComputeL2Error(exact_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 real_t 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
447real_t LogarithmGridFunctionCoefficient::Eval(ElementTransformation &T,
448 const IntegrationPoint &ip)
449{
450 MFEM_ASSERT(u != NULL, "grid function is not set");
451
452 real_t val = u->GetValue(T, ip) - obstacle->Eval(T, ip);
453 return max(min_val, log(val));
454}
455
456real_t ExponentialGridFunctionCoefficient::Eval(ElementTransformation &T,
457 const IntegrationPoint &ip)
458{
459 MFEM_ASSERT(u != NULL, "grid function is not set");
460
461 real_t val = u->GetValue(T, ip);
462 return min(max_val, max(min_val, exp(val) + obstacle->Eval(T, ip)));
463}
464
466{
467 real_t x = pt(0), y = pt(1);
468 real_t r = sqrt(x*x + y*y);
469 real_t r0 = 0.5;
470 real_t beta = 0.9;
471
472 real_t b = r0*beta;
473 real_t tmp = sqrt(r0*r0 - b*b);
474 real_t B = tmp + b*b/tmp;
475 real_t 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 real_t x = pt(0), y = pt(1);
490 real_t r = sqrt(x*x + y*y);
491 real_t r0 = 0.5;
492 real_t a = 0.348982574111686;
493 real_t 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 real_t x = pt(0), y = pt(1);
508 real_t r = sqrt(x*x + y*y);
509 real_t r0 = 0.5;
510 real_t a = 0.348982574111686;
511 real_t 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}
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
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
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.
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
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition fespace.hpp:713
A general function coefficient.
GMRES method.
Definition solvers.hpp:547
void SetKDim(int dim)
Set the number of iteration to perform between restarts, default is 50.
Definition solvers.hpp:559
virtual void Mult(const Vector &b, Vector &x) const
Iterative solution of the linear system using the GMRES method.
Definition solvers.cpp:980
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
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
The BoomerAMG solver in hypre.
Definition hypre.hpp:1691
void SetPrintLevel(int print_level)
Definition hypre.hpp:1771
Wrapper for hypre's ParCSR matrix class.
Definition hypre.hpp:388
HypreParMatrix * Transpose() const
Returns the transpose of *this.
Definition hypre.cpp:1684
Parallel smoothers in hypre.
Definition hypre.hpp:1020
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
Definition hypre.hpp:74
Class for integration point with weight.
Definition intrules.hpp:35
virtual void SetOperator(const Operator &op) override
Also calls SetOperator for the preconditioner if there is one.
Definition solvers.cpp:179
void SetRelTol(real_t rtol)
Definition solvers.hpp:209
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition solvers.cpp:173
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition solvers.cpp:71
void SetMaxIter(int max_it)
Definition solvers.hpp:211
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition fe_coll.hpp:330
Adds new Domain Integrator. Assumes ownership of lfi.
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
void Clear()
Clear the contents of the Mesh.
Definition mesh.hpp:730
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
void Assemble(int skip_zeros=1)
Adds a domain integrator. Assumes ownership of bfi.
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
static int WorldSize()
Return the size of MPI_COMM_WORLD.
static void Init(int &argc, char **&argv, int required=default_thread_required, int *provided=nullptr)
Singleton creation with Mpi::Init(argc, argv).
@ 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.
Class for parallel bilinear form.
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 Assemble(int skip_zeros=1)
Assemble the local matrix.
virtual void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)
Form the linear system matrix A, see FormLinearSystem() for details.
Abstract parallel finite element space.
Definition pfespace.hpp:29
void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1) const override
int GetTrueVSize() const override
Return the number of local vector true dofs.
Definition pfespace.hpp:289
Class for parallel grid function.
Definition pgridfunc.hpp:33
real_t ComputeH1Error(Coefficient *exsol, VectorCoefficient *exgrad, Coefficient *ell_coef, real_t Nu, int norm_type) const override
real_t ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const override
void SetFromTrueDofs(const Vector &tv) override
Set the GridFunction from the given true-dof vector.
void MakeRef(FiniteElementSpace *f, real_t *v) override
Make the ParGridFunction reference external data on a new FiniteElementSpace.
void ProjectCoefficient(Coefficient &coeff) override
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Class for parallel linear form.
void Update(ParFiniteElementSpace *pf=NULL)
Update the object according to the given new FE space *pf.
void Assemble()
Assembles the ParLinearForm i.e. sums over all domain/bdr integrators.
Class for parallel meshes.
Definition pmesh.hpp:34
Class for parallel bilinear form using different test and trial FE spaces.
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...
Scalar coefficient defined as the product of two scalar coefficients or a scalar and a scalar coeffic...
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 ex36p.cpp:505
real_t spherical_obstacle(const Vector &pt)
Definition ex36p.cpp:465
real_t exact_solution_obstacle(const Vector &pt)
Definition ex36p.cpp:487
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
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[]