MFEM  v4.3.0 Finite element discretization library
ex27p.cpp
Go to the documentation of this file.
1 // MFEM Example 27 - Parallel Version
2 //
3 // Compile with: make ex27p
4 //
5 // Sample runs: mpirun -np 4 ex27p
6 // mpirun -np 4 ex27p -dg
7 // mpirun -np 4 ex27p -dg -dbc 8 -nbc -2
8 // mpirun -np 4 ex27p -rbc-a 1 -rbc-b 8
9 //
10 // Description: This example code demonstrates the use of MFEM to define a
11 // simple finite element discretization of the Laplace problem
12 // -Delta u = 0 with a variety of boundary conditions.
13 //
14 // Specifically, we discretize using a FE space of the specified
15 // order using a continuous or discontinuous space. We then apply
16 // Dirichlet, Neumann (both homogeneous and inhomogeneous), Robin,
17 // and Periodic boundary conditions on different portions of a
18 // predefined mesh.
19 //
20 // The predefined mesh consists of a rectangle with two holes
21 // removed (see below). The narrow ends of the mesh are connected
22 // to form a Periodic boundary condition. The lower edge (tagged
23 // with attribute 1) receives an inhomogeneous Neumann boundary
24 // condition. A Robin boundary condition is applied to upper edge
25 // (attribute 2). The circular hole on the left (attribute 3)
26 // enforces a Dirichlet boundary condition. Finally, a natural
27 // boundary condition, or homogeneous Neumann BC, is applied to
28 // the circular hole on the right (attribute 4).
29 //
30 // Attribute 3 ^ y Attribute 2
31 // \ | /
32 // +-----------+-----------+
33 // | \_ | _ |
34 // | / \ | / \ |
35 // <--+---+---+---+---+---+---+--> x
36 // | \_/ | \_/ |
37 // | | \ |
38 // +-----------+-----------+ (hole radii are
39 // / | \ adjustable)
40 // Attribute 1 v Attribute 4
41 //
42 // The boundary conditions are defined as (where u is the solution
43 // field):
44 //
45 // Dirichlet: u = d
46 // Neumann: n.Grad(u) = g
47 // Robin: n.Grad(u) + a u = b
48 //
49 // The user can adjust the values of 'd', 'g', 'a', and 'b' with
50 // command line options.
51 //
52 // This example highlights the differing implementations of
53 // boundary conditions with continuous and discontinuous Galerkin
54 // formulations of the Laplace problem.
55 //
56 // We recommend viewing Examples 1 and 14 before viewing this
57 // example.
58
59 #include "mfem.hpp"
60 #include <fstream>
61 #include <iostream>
62
63 using namespace std;
64 using namespace mfem;
65
66 static double a_ = 0.2;
67
68 // Normal to hole with boundary attribute 4
69 void n4Vec(const Vector &x, Vector &n) { n = x; n[0] -= 0.5; n /= -n.Norml2(); }
70
71 Mesh * GenerateSerialMesh(int ref);
72
73 // Compute the average value of alpha*n.Grad(sol) + beta*sol over the boundary
74 // attributes marked in bdr_marker. Also computes the L2 norm of
75 // alpha*n.Grad(sol) + beta*sol - gamma over the same boundary.
76 double IntegrateBC(const ParGridFunction &sol, const Array<int> &bdr_marker,
77  double alpha, double beta, double gamma,
78  double &err);
79
80 int main(int argc, char *argv[])
81 {
82  // 1. Initialize MPI.
83  MPI_Session mpi;
84  if (!mpi.Root()) { mfem::out.Disable(); mfem::err.Disable(); }
85
86  // 2. Parse command-line options.
87  int ser_ref_levels = 2;
88  int par_ref_levels = 1;
89  int order = 1;
90  double sigma = -1.0;
91  double kappa = -1.0;
92  bool h1 = true;
93  bool visualization = true;
94
95  double mat_val = 1.0;
96  double dbc_val = 0.0;
97  double nbc_val = 1.0;
98  double rbc_a_val = 1.0; // du/dn + a * u = b
99  double rbc_b_val = 1.0;
100
101  OptionsParser args(argc, argv);
102  args.AddOption(&h1, "-h1", "--continuous", "-dg", "--discontinuous",
103  "Select continuous \"H1\" or discontinuous \"DG\" basis.");
104  args.AddOption(&order, "-o", "--order",
105  "Finite element order (polynomial degree) or -1 for"
106  " isoparametric space.");
107  args.AddOption(&sigma, "-s", "--sigma",
108  "One of the two DG penalty parameters, typically +1/-1."
109  " See the documentation of class DGDiffusionIntegrator.");
110  args.AddOption(&kappa, "-k", "--kappa",
111  "One of the two DG penalty parameters, should be positive."
112  " Negative values are replaced with (order+1)^2.");
113  args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
114  "Number of times to refine the mesh uniformly in serial.");
115  args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
116  "Number of times to refine the mesh uniformly in parallel.");
117  args.AddOption(&mat_val, "-mat", "--material-value",
118  "Constant value for material coefficient "
119  "in the Laplace operator.");
120  args.AddOption(&dbc_val, "-dbc", "--dirichlet-value",
121  "Constant value for Dirichlet Boundary Condition.");
122  args.AddOption(&nbc_val, "-nbc", "--neumann-value",
123  "Constant value for Neumann Boundary Condition.");
124  args.AddOption(&rbc_a_val, "-rbc-a", "--robin-a-value",
125  "Constant 'a' value for Robin Boundary Condition: "
126  "du/dn + a * u = b.");
127  args.AddOption(&rbc_b_val, "-rbc-b", "--robin-b-value",
128  "Constant 'b' value for Robin Boundary Condition: "
129  "du/dn + a * u = b.");
131  "Radius of holes in the mesh.");
132  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
133  "--no-visualization",
134  "Enable or disable GLVis visualization.");
135  args.Parse();
136  if (!args.Good())
137  {
138  args.PrintUsage(mfem::out);
139  return 1;
140  }
141  if (kappa < 0 && !h1)
142  {
143  kappa = (order+1)*(order+1);
144  }
145  args.PrintOptions(mfem::out);
146
147  if (a_ < 0.01)
148  {
149  mfem::out << "Hole radius too small, resetting to 0.01.\n";
150  a_ = 0.01;
151  }
152  if (a_ > 0.49)
153  {
154  mfem::out << "Hole radius too large, resetting to 0.49.\n";
155  a_ = 0.49;
156  }
157
158  // 3. Construct the (serial) mesh and refine it if requested.
159  Mesh *mesh = GenerateSerialMesh(ser_ref_levels);
160  int dim = mesh->Dimension();
161
162  // 4. Define a parallel mesh by a partitioning of the serial mesh. Refine
163  // this mesh further in parallel to increase the resolution. Once the
164  // parallel mesh is defined, the serial mesh can be deleted.
165  ParMesh pmesh(MPI_COMM_WORLD, *mesh);
166  delete mesh;
167  for (int l = 0; l < par_ref_levels; l++)
168  {
169  pmesh.UniformRefinement();
170  }
171
172  // 5. Define a parallel finite element space on the parallel mesh. Here we
173  // use either continuous Lagrange finite elements or discontinuous
174  // Galerkin finite elements of the specified order.
176  h1 ? (FiniteElementCollection*)new H1_FECollection(order, dim) :
177  (FiniteElementCollection*)new DG_FECollection(order, dim);
178  ParFiniteElementSpace fespace(&pmesh, fec);
179  HYPRE_BigInt size = fespace.GlobalTrueVSize();
180  mfem::out << "Number of finite element unknowns: " << size << endl;
181
182  // 6. Create "marker arrays" to define the portions of boundary associated
183  // with each type of boundary condition. These arrays have an entry
184  // corresponding to each boundary attribute. Placing a '1' in entry i
185  // marks attribute i+1 as being active, '0' is inactive.
186  Array<int> nbc_bdr(pmesh.bdr_attributes.Max());
187  Array<int> rbc_bdr(pmesh.bdr_attributes.Max());
188  Array<int> dbc_bdr(pmesh.bdr_attributes.Max());
189
190  nbc_bdr = 0; nbc_bdr[0] = 1;
191  rbc_bdr = 0; rbc_bdr[1] = 1;
192  dbc_bdr = 0; dbc_bdr[2] = 1;
193
194  Array<int> ess_tdof_list(0);
195  if (h1 && pmesh.bdr_attributes.Size())
196  {
197  // For a continuous basis the linear system must be modifed to enforce an
198  // essential (Dirichlet) boundary condition. In the DG case this is not
199  // necessary as the boundary condition will only be enforced weakly.
200  fespace.GetEssentialTrueDofs(dbc_bdr, ess_tdof_list);
201  }
202
203  // 7. Setup the various coefficients needed for the Laplace operator and the
204  // various boundary conditions. In general these coefficients could be
205  // functions of position but here we use only constants.
206  ConstantCoefficient matCoef(mat_val);
207  ConstantCoefficient dbcCoef(dbc_val);
208  ConstantCoefficient nbcCoef(nbc_val);
209  ConstantCoefficient rbcACoef(rbc_a_val);
210  ConstantCoefficient rbcBCoef(rbc_b_val);
211
212  // Since the n.Grad(u) terms arise by integrating -Div(m Grad(u)) by parts we
213  // must introduce the coefficient 'm' into the boundary conditions.
214  // Therefore, in the case of the Neumann BC, we actually enforce m n.Grad(u)
215  // = m g rather than simply n.Grad(u) = g.
216  ProductCoefficient m_nbcCoef(matCoef, nbcCoef);
217  ProductCoefficient m_rbcACoef(matCoef, rbcACoef);
218  ProductCoefficient m_rbcBCoef(matCoef, rbcBCoef);
219
220  // 8. Define the solution vector u as a parallel finite element grid function
221  // corresponding to fespace. Initialize u with initial guess of zero.
222  ParGridFunction u(&fespace);
223  u = 0.0;
224
225  // 9. Set up the parallel bilinear form a(.,.) on the finite element space
226  // corresponding to the Laplacian operator -Delta, by adding the Diffusion
227  // domain integrator.
228  ParBilinearForm a(&fespace);
230  if (h1)
231  {
232  // Add a Mass integrator on the Robin boundary
233  a.AddBoundaryIntegrator(new MassIntegrator(m_rbcACoef), rbc_bdr);
234  }
235  else
236  {
237  // Add the interfacial portion of the Laplace operator
239  sigma, kappa));
240
241  // Counteract the n.Grad(u) term on the Dirichlet portion of the boundary
242  a.AddBdrFaceIntegrator(new DGDiffusionIntegrator(matCoef, sigma, kappa),
243  dbc_bdr);
244
245  // Augment the n.Grad(u) term with a*u on the Robin portion of boundary
247  rbc_bdr);
248  }
249  a.Assemble();
250
251  // 10. Assemble the parallel linear form for the right hand side vector.
252  ParLinearForm b(&fespace);
253
254  if (h1)
255  {
256  // Set the Dirichlet values in the solution vector
257  u.ProjectBdrCoefficient(dbcCoef, dbc_bdr);
258
259  // Add the desired value for n.Grad(u) on the Neumann boundary
260  b.AddBoundaryIntegrator(new BoundaryLFIntegrator(m_nbcCoef), nbc_bdr);
261
262  // Add the desired value for n.Grad(u) + a*u on the Robin boundary
263  b.AddBoundaryIntegrator(new BoundaryLFIntegrator(m_rbcBCoef), rbc_bdr);
264  }
265  else
266  {
267  // Add the desired value for the Dirichlet boundary
268  b.AddBdrFaceIntegrator(new DGDirichletLFIntegrator(dbcCoef, matCoef,
269  sigma, kappa),
270  dbc_bdr);
271
272  // Add the desired value for n.Grad(u) on the Neumann boundary
274  nbc_bdr);
275
276  // Add the desired value for n.Grad(u) + a*u on the Robin boundary
278  rbc_bdr);
279  }
280  b.Assemble();
281
282  // 11. Construct the linear system.
283  OperatorPtr A;
284  Vector B, X;
285  a.FormLinearSystem(ess_tdof_list, u, b, A, X, B);
286
287  // 12. Solve the linear system A X = B.
288  HypreSolver *amg = new HypreBoomerAMG;
289  if (h1 || sigma == -1.0)
290  {
291  HyprePCG pcg(MPI_COMM_WORLD);
292  pcg.SetTol(1e-12);
293  pcg.SetMaxIter(200);
294  pcg.SetPrintLevel(2);
295  pcg.SetPreconditioner(*amg);
296  pcg.SetOperator(*A);
297  pcg.Mult(B, X);
298  }
299  else
300  {
301  GMRESSolver gmres(MPI_COMM_WORLD);
302  gmres.SetAbsTol(0.0);
303  gmres.SetRelTol(1e-12);
304  gmres.SetMaxIter(200);
305  gmres.SetKDim(10);
306  gmres.SetPrintLevel(1);
307  gmres.SetPreconditioner(*amg);
308  gmres.SetOperator(*A);
309  gmres.Mult(B, X);
310  }
311  delete amg;
312
313  // 13. Recover the parallel grid function corresponding to U. This is the
314  // local finite element solution on each processor.
315  a.RecoverFEMSolution(X, b, u);
316
317  // 14. Build a mass matrix to help solve for n.Grad(u) where 'n' is a surface
318  // normal.
319  ParBilinearForm m(&fespace);
321  m.Assemble();
322
323  ess_tdof_list.SetSize(0);
324  OperatorPtr M;
325  m.FormSystemMatrix(ess_tdof_list, M);
326
327  // 15. Compute the various boundary integrals.
328  mfem::out << endl
329  << "Verifying boundary conditions" << endl
330  << "=============================" << endl;
331  {
332  // Integrate the solution on the Dirichlet boundary and compare to the
333  // expected value.
334  double err, avg = IntegrateBC(u, dbc_bdr, 0.0, 1.0, dbc_val, err);
335
336  bool hom_dbc = (dbc_val == 0.0);
337  err /= hom_dbc ? 1.0 : fabs(dbc_val);
338  mfem::out << "Average of solution on Gamma_dbc:\t"
339  << avg << ", \t"
340  << (hom_dbc ? "absolute" : "relative")
341  << " error " << err << endl;
342  }
343  {
344  // Integrate n.Grad(u) on the inhomogeneous Neumann boundary and compare
345  // to the expected value.
346  double err, avg = IntegrateBC(u, nbc_bdr, 1.0, 0.0, nbc_val, err);
347
348  bool hom_nbc = (nbc_val == 0.0);
349  err /= hom_nbc ? 1.0 : fabs(nbc_val);
350  mfem::out << "Average of n.Grad(u) on Gamma_nbc:\t"
351  << avg << ", \t"
352  << (hom_nbc ? "absolute" : "relative")
353  << " error " << err << endl;
354  }
355  {
356  // Integrate n.Grad(u) on the homogeneous Neumann boundary and compare to
357  // the expected value of zero.
358  Array<int> nbc0_bdr(pmesh.bdr_attributes.Max());
359  nbc0_bdr = 0;
360  nbc0_bdr[3] = 1;
361
362  double err, avg = IntegrateBC(u, nbc0_bdr, 1.0, 0.0, 0.0, err);
363
364  bool hom_nbc = true;
365  mfem::out << "Average of n.Grad(u) on Gamma_nbc0:\t"
366  << avg << ", \t"
367  << (hom_nbc ? "absolute" : "relative")
368  << " error " << err << endl;
369  }
370  {
371  // Integrate n.Grad(u) + a * u on the Robin boundary and compare to the
372  // expected value.
373  double err, avg = IntegrateBC(u, rbc_bdr, 1.0, rbc_a_val, rbc_b_val, err);
374
375  bool hom_rbc = (rbc_b_val == 0.0);
376  err /= hom_rbc ? 1.0 : fabs(rbc_b_val);
377  mfem::out << "Average of n.Grad(u)+a*u on Gamma_rbc:\t"
378  << avg << ", \t"
379  << (hom_rbc ? "absolute" : "relative")
380  << " error " << err << endl;
381  }
382
383  // 16. Save the refined mesh and the solution in parallel. This output can be
384  // viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
385  {
386  ostringstream mesh_name, sol_name;
387  mesh_name << "mesh." << setfill('0') << setw(6) << mpi.WorldRank();
388  sol_name << "sol." << setfill('0') << setw(6) << mpi.WorldRank();
389
390  ofstream mesh_ofs(mesh_name.str().c_str());
391  mesh_ofs.precision(8);
392  pmesh.Print(mesh_ofs);
393
394  ofstream sol_ofs(sol_name.str().c_str());
395  sol_ofs.precision(8);
396  u.Save(sol_ofs);
397  }
398
399  // 17. Send the solution by socket to a GLVis server.
400  if (visualization)
401  {
402  string title_str = h1 ? "H1" : "DG";
403  char vishost[] = "localhost";
404  int visport = 19916;
405  socketstream sol_sock(vishost, visport);
406  sol_sock << "parallel " << mpi.WorldSize()
407  << " " << mpi.WorldRank() << "\n";
408  sol_sock.precision(8);
409  sol_sock << "solution\n" << pmesh << u
410  << "window_title '" << title_str << " Solution'"
411  << " keys 'mmc'" << flush;
412  }
413
414  // 18. Free the used memory.
415  delete fec;
416
417  return 0;
418 }
419
420 void quad_trans(double u, double v, double &x, double &y, bool log = false)
421 {
422  double a = a_; // Radius of disc
423
424  double d = 4.0 * a * (M_SQRT2 - 2.0 * a) * (1.0 - 2.0 * v);
425
426  double v0 = (1.0 + M_SQRT2) * (M_SQRT2 * a - 2.0 * v) *
427  ((4.0 - 3 * M_SQRT2) * a +
428  (8.0 * (M_SQRT2 - 1.0) * a - 2.0) * v) / d;
429
430  double r = 2.0 * ((M_SQRT2 - 1.0) * a * a * (1.0 - 4.0 *v) +
431  2.0 * (1.0 + M_SQRT2 *
432  (1.0 + 2.0 * (2.0 * a - M_SQRT2 - 1.0) * a)) * v * v
433  ) / d;
434
435  double t = asin(v / r) * u / v;
436  if (log)
437  {
438  mfem::out << "u, v, r, v0, t "
439  << u << " " << v << " " << r << " " << v0 << " " << t
440  << endl;
441  }
442  x = r * sin(t);
443  y = r * cos(t) - v0;
444 }
445
446 void trans(const Vector &u, Vector &x)
447 {
448  double tol = 1e-4;
449
450  if (u[1] > 0.5 - tol || u[1] < -0.5 + tol)
451  {
452  x = u;
453  return;
454  }
455  if (u[0] > 1.0 - tol || u[0] < -1.0 + tol || fabs(u[0]) < tol)
456  {
457  x = u;
458  return;
459  }
460
461  if (u[0] > 0.0)
462  {
463  if (u[1] > fabs(u[0] - 0.5))
464  {
465  quad_trans(u[0] - 0.5, u[1], x[0], x[1]);
466  x[0] += 0.5;
467  return;
468  }
469  if (u[1] < -fabs(u[0] - 0.5))
470  {
471  quad_trans(u[0] - 0.5, -u[1], x[0], x[1]);
472  x[0] += 0.5;
473  x[1] *= -1.0;
474  return;
475  }
476  if (u[0] - 0.5 > fabs(u[1]))
477  {
478  quad_trans(u[1], u[0] - 0.5, x[1], x[0]);
479  x[0] += 0.5;
480  return;
481  }
482  if (u[0] - 0.5 < -fabs(u[1]))
483  {
484  quad_trans(u[1], 0.5 - u[0], x[1], x[0]);
485  x[0] *= -1.0;
486  x[0] += 0.5;
487  return;
488  }
489  }
490  else
491  {
492  if (u[1] > fabs(u[0] + 0.5))
493  {
494  quad_trans(u[0] + 0.5, u[1], x[0], x[1]);
495  x[0] -= 0.5;
496  return;
497  }
498  if (u[1] < -fabs(u[0] + 0.5))
499  {
500  quad_trans(u[0] + 0.5, -u[1], x[0], x[1]);
501  x[0] -= 0.5;
502  x[1] *= -1.0;
503  return;
504  }
505  if (u[0] + 0.5 > fabs(u[1]))
506  {
507  quad_trans(u[1], u[0] + 0.5, x[1], x[0]);
508  x[0] -= 0.5;
509  return;
510  }
511  if (u[0] + 0.5 < -fabs(u[1]))
512  {
513  quad_trans(u[1], -0.5 - u[0], x[1], x[0]);
514  x[0] *= -1.0;
515  x[0] -= 0.5;
516  return;
517  }
518  }
519  x = u;
520 }
521
522 Mesh * GenerateSerialMesh(int ref)
523 {
524  Mesh * mesh = new Mesh(2, 29, 16, 24, 2);
525
526  int vi[4];
527
528  for (int i=0; i<2; i++)
529  {
530  int o = 13 * i;
531  vi[0] = o + 0; vi[1] = o + 3; vi[2] = o + 4; vi[3] = o + 1;
533
534  vi[0] = o + 1; vi[1] = o + 4; vi[2] = o + 5; vi[3] = o + 2;
536
537  vi[0] = o + 5; vi[1] = o + 8; vi[2] = o + 9; vi[3] = o + 2;
539
540  vi[0] = o + 8; vi[1] = o + 12; vi[2] = o + 15; vi[3] = o + 9;
542
543  vi[0] = o + 11; vi[1] = o + 14; vi[2] = o + 15; vi[3] = o + 12;
545
546  vi[0] = o + 10; vi[1] = o + 13; vi[2] = o + 14; vi[3] = o + 11;
548
549  vi[0] = o + 6; vi[1] = o + 13; vi[2] = o + 10; vi[3] = o + 7;
551
552  vi[0] = o + 0; vi[1] = o + 6; vi[2] = o + 7; vi[3] = o + 3;
554  }
555
556  vi[0] = 0; vi[1] = 6; mesh->AddBdrSegment(vi, 1);
557  vi[0] = 6; vi[1] = 13; mesh->AddBdrSegment(vi, 1);
558  vi[0] = 13; vi[1] = 19; mesh->AddBdrSegment(vi, 1);
559  vi[0] = 19; vi[1] = 26; mesh->AddBdrSegment(vi, 1);
560
561  vi[0] = 28; vi[1] = 22; mesh->AddBdrSegment(vi, 2);
562  vi[0] = 22; vi[1] = 15; mesh->AddBdrSegment(vi, 2);
563  vi[0] = 15; vi[1] = 9; mesh->AddBdrSegment(vi, 2);
564  vi[0] = 9; vi[1] = 2; mesh->AddBdrSegment(vi, 2);
565
566  for (int i=0; i<2; i++)
567  {
568  int o = 13 * i;
569  vi[0] = o + 7; vi[1] = o + 3; mesh->AddBdrSegment(vi, 3 + i);
570  vi[0] = o + 10; vi[1] = o + 7; mesh->AddBdrSegment(vi, 3 + i);
571  vi[0] = o + 11; vi[1] = o + 10; mesh->AddBdrSegment(vi, 3 + i);
572  vi[0] = o + 12; vi[1] = o + 11; mesh->AddBdrSegment(vi, 3 + i);
573  vi[0] = o + 8; vi[1] = o + 12; mesh->AddBdrSegment(vi, 3 + i);
574  vi[0] = o + 5; vi[1] = o + 8; mesh->AddBdrSegment(vi, 3 + i);
575  vi[0] = o + 4; vi[1] = o + 5; mesh->AddBdrSegment(vi, 3 + i);
576  vi[0] = o + 3; vi[1] = o + 4; mesh->AddBdrSegment(vi, 3 + i);
577  }
578
579  double d[2];
580  double a = a_ / M_SQRT2;
581
582  d[0] = -1.0; d[1] = -0.5; mesh->AddVertex(d);
583  d[0] = -1.0; d[1] = 0.0; mesh->AddVertex(d);
584  d[0] = -1.0; d[1] = 0.5; mesh->AddVertex(d);
585
586  d[0] = -0.5 - a; d[1] = -a; mesh->AddVertex(d);
587  d[0] = -0.5 - a; d[1] = 0.0; mesh->AddVertex(d);
588  d[0] = -0.5 - a; d[1] = a; mesh->AddVertex(d);
589
590  d[0] = -0.5; d[1] = -0.5; mesh->AddVertex(d);
591  d[0] = -0.5; d[1] = -a; mesh->AddVertex(d);
592  d[0] = -0.5; d[1] = a; mesh->AddVertex(d);
593  d[0] = -0.5; d[1] = 0.5; mesh->AddVertex(d);
594
595  d[0] = -0.5 + a; d[1] = -a; mesh->AddVertex(d);
596  d[0] = -0.5 + a; d[1] = 0.0; mesh->AddVertex(d);
597  d[0] = -0.5 + a; d[1] = a; mesh->AddVertex(d);
598
599  d[0] = 0.0; d[1] = -0.5; mesh->AddVertex(d);
600  d[0] = 0.0; d[1] = 0.0; mesh->AddVertex(d);
601  d[0] = 0.0; d[1] = 0.5; mesh->AddVertex(d);
602
603  d[0] = 0.5 - a; d[1] = -a; mesh->AddVertex(d);
604  d[0] = 0.5 - a; d[1] = 0.0; mesh->AddVertex(d);
605  d[0] = 0.5 - a; d[1] = a; mesh->AddVertex(d);
606
607  d[0] = 0.5; d[1] = -0.5; mesh->AddVertex(d);
608  d[0] = 0.5; d[1] = -a; mesh->AddVertex(d);
609  d[0] = 0.5; d[1] = a; mesh->AddVertex(d);
610  d[0] = 0.5; d[1] = 0.5; mesh->AddVertex(d);
611
612  d[0] = 0.5 + a; d[1] = -a; mesh->AddVertex(d);
613  d[0] = 0.5 + a; d[1] = 0.0; mesh->AddVertex(d);
614  d[0] = 0.5 + a; d[1] = a; mesh->AddVertex(d);
615
616  d[0] = 1.0; d[1] = -0.5; mesh->AddVertex(d);
617  d[0] = 1.0; d[1] = 0.0; mesh->AddVertex(d);
618  d[0] = 1.0; d[1] = 0.5; mesh->AddVertex(d);
619
620  mesh->FinalizeTopology();
621
622  mesh->SetCurvature(1, true);
623
624  // Stitch the ends of the stack together
625  {
626  Array<int> v2v(mesh->GetNV());
627  for (int i = 0; i < v2v.Size() - 3; i++)
628  {
629  v2v[i] = i;
630  }
631  // identify vertices on the narrow ends of the rectangle
632  v2v[v2v.Size() - 3] = 0;
633  v2v[v2v.Size() - 2] = 1;
634  v2v[v2v.Size() - 1] = 2;
635
636  // renumber elements
637  for (int i = 0; i < mesh->GetNE(); i++)
638  {
639  Element *el = mesh->GetElement(i);
640  int *v = el->GetVertices();
641  int nv = el->GetNVertices();
642  for (int j = 0; j < nv; j++)
643  {
644  v[j] = v2v[v[j]];
645  }
646  }
647  // renumber boundary elements
648  for (int i = 0; i < mesh->GetNBE(); i++)
649  {
650  Element *el = mesh->GetBdrElement(i);
651  int *v = el->GetVertices();
652  int nv = el->GetNVertices();
653  for (int j = 0; j < nv; j++)
654  {
655  v[j] = v2v[v[j]];
656  }
657  }
658  mesh->RemoveUnusedVertices();
659  mesh->RemoveInternalBoundaries();
660  }
661  mesh->SetCurvature(3, true);
662
663  for (int l = 0; l < ref; l++)
664  {
665  mesh->UniformRefinement();
666  }
667
668  mesh->Transform(trans);
669
670  return mesh;
671 }
672
673 double IntegrateBC(const ParGridFunction &x, const Array<int> &bdr,
674  double alpha, double beta, double gamma,
675  double &glb_err)
676 {
677  double loc_vals[3];
678  double &nrm = loc_vals[0];
679  double &avg = loc_vals[1];
680  double &err = loc_vals[2];
681
682  nrm = 0.0;
683  avg = 0.0;
684  err = 0.0;
685
686  const bool a_is_zero = alpha == 0.0;
687  const bool b_is_zero = beta == 0.0;
688
689  const ParFiniteElementSpace &fes = *x.ParFESpace();
690  MFEM_ASSERT(fes.GetVDim() == 1, "");
691  ParMesh &mesh = *fes.GetParMesh();
692  Vector shape, loc_dofs, w_nor;
693  DenseMatrix dshape;
694  Array<int> dof_ids;
695  for (int i = 0; i < mesh.GetNBE(); i++)
696  {
697  if (bdr[mesh.GetBdrAttribute(i)-1] == 0) { continue; }
698
699  FaceElementTransformations *FTr = mesh.GetBdrFaceTransformations(i);
700  if (FTr == nullptr) { continue; }
701
702  const FiniteElement &fe = *fes.GetFE(FTr->Elem1No);
703  MFEM_ASSERT(fe.GetMapType() == FiniteElement::VALUE, "");
704  const int int_order = 2*fe.GetOrder() + 3;
705  const IntegrationRule &ir = IntRules.Get(FTr->FaceGeom, int_order);
706
707  fes.GetElementDofs(FTr->Elem1No, dof_ids);
708  x.GetSubVector(dof_ids, loc_dofs);
709  if (!a_is_zero)
710  {
711  const int sdim = FTr->Face->GetSpaceDim();
712  w_nor.SetSize(sdim);
713  dshape.SetSize(fe.GetDof(), sdim);
714  }
715  if (!b_is_zero)
716  {
717  shape.SetSize(fe.GetDof());
718  }
719  for (int j = 0; j < ir.GetNPoints(); j++)
720  {
721  const IntegrationPoint &ip = ir.IntPoint(j);
722  IntegrationPoint eip;
723  FTr->Loc1.Transform(ip, eip);
724  FTr->Face->SetIntPoint(&ip);
725  double face_weight = FTr->Face->Weight();
726  double val = 0.0;
727  if (!a_is_zero)
728  {
729  FTr->Elem1->SetIntPoint(&eip);
730  fe.CalcPhysDShape(*FTr->Elem1, dshape);
731  CalcOrtho(FTr->Face->Jacobian(), w_nor);
732  val += alpha * dshape.InnerProduct(w_nor, loc_dofs) / face_weight;
733  }
734  if (!b_is_zero)
735  {
736  fe.CalcShape(eip, shape);
737  val += beta * (shape * loc_dofs);
738  }
739
740  // Measure the length of the boundary
741  nrm += ip.weight * face_weight;
742
743  // Integrate alpha * n.Grad(x) + beta * x
744  avg += val * ip.weight * face_weight;
745
746  // Integrate |alpha * n.Grad(x) + beta * x - gamma|^2
747  val -= gamma;
748  err += (val*val) * ip.weight * face_weight;
749  }
750  }
751
752  double glb_vals[3];
753  MPI_Allreduce(loc_vals, glb_vals, 3, MPI_DOUBLE, MPI_SUM, fes.GetComm());
754
755  double glb_nrm = glb_vals[0];
756  double glb_avg = glb_vals[1];
757  glb_err = glb_vals[2];
758
759  // Normalize by the length of the boundary
760  if (std::abs(glb_nrm) > 0.0)
761  {
762  glb_err /= glb_nrm;
763  glb_avg /= glb_nrm;
764  }
765
766  // Compute l2 norm of the error in the boundary condition (negative
767  // quadrature weights may produce negative 'err')
768  glb_err = (glb_err >= 0.0) ? sqrt(glb_err) : -sqrt(-glb_err);
769
770  // Return the average value of alpha * n.Grad(x) + beta * x
771  return glb_avg;
772 }
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247
Abstract class for all finite elements.
Definition: fe.hpp:243
int WorldSize() const
Return MPI_COMM_WORLD&#39;s size.
void SetTol(double tol)
Definition: hypre.cpp:3569
int Size() const
Return the logical size of the array.
Definition: array.hpp:134
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Definition: pfespace.cpp:790
void trans(const Vector &u, Vector &x)
Definition: ex27.cpp:421
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
int AddQuad(int v1, int v2, int v3, int v4, int attr=1)
Definition: mesh.cpp:1324
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:920
ParMesh * GetParMesh() const
Definition: pfespace.hpp:267
A coefficient that is constant across space and time.
Definition: coefficient.hpp:78
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(...
virtual void GetVertices(Array< int > &v) const =0
Returns element&#39;s vertices.
void n4Vec(const Vector &x, Vector &n)
Definition: ex27.cpp:69
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: hypre.cpp:3547
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:513
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:849
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:783
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
Scalar coefficient defined as the product of two scalar coefficients or a scalar and a scalar coeffic...
double InnerProduct(const double *x, const double *y) const
Compute y^t A x.
Definition: densemat.cpp:299
void quad_trans(double u, double v, double &x, double &y, bool log=false)
Definition: ex27.cpp:395
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition: eltrans.hpp:85
HYPRE_BigInt GlobalTrueVSize() const
Definition: pfespace.hpp:275
A specialized ElementTransformation class representing a face and its two neighboring elements...
Definition: eltrans.hpp:467
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition: vector.cpp:525
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe.hpp:327
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
Class for boundary integration L(v) := (g, v)
Definition: lininteg.hpp:154
void Transform(void(*f)(const Vector &, Vector &))
Definition: mesh.cpp:11008
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:846
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:841
double kappa
Definition: ex24.cpp:54
Abstract parallel finite element space.
Definition: pfespace.hpp:28
int GetMapType() const
Returns the FiniteElement::MapType of the element describing how reference functions are mapped to ph...
Definition: fe.hpp:349
void CalcOrtho(const DenseMatrix &J, Vector &n)
Definition: densemat.cpp:2303
void RemoveInternalBoundaries()
Definition: mesh.cpp:11164
void SetPrintLevel(int print_lvl)
Definition: hypre.cpp:3589
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:884
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1387
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:250
int AddVertex(double x, double y=0.0, double z=0.0)
Definition: mesh.cpp:1263
Class for parallel linear form.
Definition: plinearform.hpp:26
A simple convenience class that calls MPI_Init() at construction and MPI_Finalize() at destruction...
void SetPrintLevel(int print_lvl)
Definition: solvers.cpp:71
MPI_Comm GetComm() const
Definition: pfespace.hpp:263
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:98
virtual const FiniteElement * GetFE(int i) const
Definition: pfespace.cpp:511
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:150
constexpr char vishost[]
IntegrationPointTransformation Loc1
Definition: eltrans.hpp:511
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9143
constexpr int visport
const DenseMatrix & Jacobian()
Return the Jacobian matrix of the transformation at the currently set IntegrationPoint, using the method SetIntPoint().
Definition: eltrans.hpp:111
void SetMaxIter(int max_it)
Definition: solvers.hpp:100
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:4882
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const =0
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
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:361
bool Root() const
Return true if WorldRank() == 0.
virtual void Print(std::ostream &out=mfem::out) const
Definition: pmesh.cpp:4382
int AddBdrSegment(int v1, int v2, int attr=1)
Definition: mesh.cpp:1440
double Weight()
Return the weight of the Jacobian matrix of the transformation at the currently set IntegrationPoint...
Definition: eltrans.hpp:123
const Element * GetElement(int i) const
Definition: mesh.hpp:942
Adds new Boundary Face Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:83
void CalcPhysDShape(ElementTransformation &Trans, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in physical space at the poi...
Definition: fe.cpp:192
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:534
Adds new Boundary Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:70
int Dimension() const
Definition: mesh.hpp:911
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:457
void SetMaxIter(int max_iter)
Definition: hypre.cpp:3579
L2_FECollection DG_FECollection
Declare an alternative name for L2_FECollection = DG_FECollection.
Definition: fe_coll.hpp:338
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
int WorldRank() const
Return MPI_COMM_WORLD&#39;s rank.
void SetRelTol(double rtol)
Definition: solvers.hpp:98
PCG solver in hypre.
Definition: hypre.hpp:1060
void ProjectBdrCoefficient(Coefficient *coeff[], VectorCoefficient *vcoeff, Array< int > &attr)
Definition: pgridfunc.cpp:601
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe.hpp:323
OutStream err(std::cerr)
Global stream used by the library for standard error output. Initially it uses the same std::streambu...
Definition: globals.hpp:71
void RemoveUnusedVertices()
Remove unused vertices and rebuild mesh connectivity.
Definition: mesh.cpp:11057
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:26
void FinalizeTopology(bool generate_bdr=true)
Finalize the construction of the secondary topology (connectivity) data of a Mesh.
Definition: mesh.cpp:2507
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:674
HYPRE_Int HYPRE_BigInt
void Transform(const IntegrationPoint &, IntegrationPoint &)
Definition: eltrans.cpp:533
GMRES method.
Definition: solvers.hpp:348
void Disable()
Disable output.
Definition: globals.hpp:54
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
virtual int GetSpaceDim() const =0
Get the dimension of the target (physical) space.
double a
Definition: lissajous.cpp:41
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
Definition: mesh.hpp:843
virtual void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)
Form the linear system matrix A, see FormLinearSystem() for details.
virtual void GetElementDofs(int i, Array< int > &dofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: pfespace.cpp:467
Class for integration point with weight.
Definition: intrules.hpp:25
Adds new Boundary Integrator. Assumes ownership of bfi.
int dim
Definition: ex24.cpp:53
Adds new interior Face Integrator. Assumes ownership of bfi.
void SetPreconditioner(HypreSolver &precond)
Set the hypre solver to be used as a preconditioner.
Definition: hypre.cpp:3594
Adds new Domain Integrator. Assumes ownership of bfi.
Mesh * GenerateSerialMesh(int ref)
Definition: ex27.cpp:497
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:327
ElementTransformation * Elem1
Definition: eltrans.hpp:509
ElementTransformation * Face
Definition: eltrans.hpp:510
virtual int GetNVertices() const =0
Class for parallel bilinear form.
Abstract class for hypre&#39;s solvers and preconditioners.
Definition: hypre.hpp:970
RefCoord t[3]
const double alpha
Definition: ex15.cpp:369
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
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
Class for parallel grid function.
Definition: pgridfunc.hpp:32
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
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
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Solve Ax=b with hypre&#39;s PCG.
Definition: hypre.cpp:3617
Class for parallel meshes.
Definition: pmesh.hpp:32
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:377
double IntegrateBC(const GridFunction &sol, const Array< int > &bdr_marker, double alpha, double beta, double gamma, double &err)
Definition: ex27.cpp:648
Abstract data type element.
Definition: element.hpp:28
int main()
const Element * GetBdrElement(int i) const
Definition: mesh.hpp:946