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