MFEM  v4.4.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
ex27.cpp
Go to the documentation of this file.
1 // MFEM Example 27 - Serial Version
2 //
3 // Compile with: make ex27
4 //
5 // Sample runs: ex27
6 // ex27 -dg
7 // ex27 -dg -dbc 8 -nbc -2
8 // ex27 -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 GridFunction &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. Parse command-line options.
83  int ser_ref_levels = 2;
84  int order = 1;
85  double sigma = -1.0;
86  double kappa = -1.0;
87  bool h1 = true;
88  bool visualization = true;
89 
90  double mat_val = 1.0;
91  double dbc_val = 0.0;
92  double nbc_val = 1.0;
93  double rbc_a_val = 1.0; // du/dn + a * u = b
94  double rbc_b_val = 1.0;
95 
96  OptionsParser args(argc, argv);
97  args.AddOption(&h1, "-h1", "--continuous", "-dg", "--discontinuous",
98  "Select continuous \"H1\" or discontinuous \"DG\" basis.");
99  args.AddOption(&order, "-o", "--order",
100  "Finite element order (polynomial degree) or -1 for"
101  " isoparametric space.");
102  args.AddOption(&sigma, "-s", "--sigma",
103  "One of the two DG penalty parameters, typically +1/-1."
104  " See the documentation of class DGDiffusionIntegrator.");
105  args.AddOption(&kappa, "-k", "--kappa",
106  "One of the two DG penalty parameters, should be positive."
107  " Negative values are replaced with (order+1)^2.");
108  args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
109  "Number of times to refine the mesh uniformly in serial.");
110  args.AddOption(&mat_val, "-mat", "--material-value",
111  "Constant value for material coefficient "
112  "in the Laplace operator.");
113  args.AddOption(&dbc_val, "-dbc", "--dirichlet-value",
114  "Constant value for Dirichlet Boundary Condition.");
115  args.AddOption(&nbc_val, "-nbc", "--neumann-value",
116  "Constant value for Neumann Boundary Condition.");
117  args.AddOption(&rbc_a_val, "-rbc-a", "--robin-a-value",
118  "Constant 'a' value for Robin Boundary Condition: "
119  "du/dn + a * u = b.");
120  args.AddOption(&rbc_b_val, "-rbc-b", "--robin-b-value",
121  "Constant 'b' value for Robin Boundary Condition: "
122  "du/dn + a * u = b.");
123  args.AddOption(&a_, "-a", "--radius",
124  "Radius of holes in the mesh.");
125  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
126  "--no-visualization",
127  "Enable or disable GLVis visualization.");
128  args.Parse();
129  if (!args.Good())
130  {
131  args.PrintUsage(mfem::out);
132  return 1;
133  }
134  if (kappa < 0 && !h1)
135  {
136  kappa = (order+1)*(order+1);
137  }
138  args.PrintOptions(mfem::out);
139 
140  if (a_ < 0.01)
141  {
142  mfem::out << "Hole radius too small, resetting to 0.01.\n";
143  a_ = 0.01;
144  }
145  if (a_ > 0.49)
146  {
147  mfem::out << "Hole radius too large, resetting to 0.49.\n";
148  a_ = 0.49;
149  }
150 
151  // 2. Construct the (serial) mesh and refine it if requested.
152  Mesh *mesh = GenerateSerialMesh(ser_ref_levels);
153  int dim = mesh->Dimension();
154 
155  // 3. Define a finite element space on the serial mesh. Here we use either
156  // continuous Lagrange finite elements or discontinuous Galerkin finite
157  // elements of the specified order.
159  h1 ? (FiniteElementCollection*)new H1_FECollection(order, dim) :
160  (FiniteElementCollection*)new DG_FECollection(order, dim);
161  FiniteElementSpace fespace(mesh, fec);
162  int size = fespace.GetTrueVSize();
163  mfem::out << "Number of finite element unknowns: " << size << endl;
164 
165  // 4. Create "marker arrays" to define the portions of boundary associated
166  // with each type of boundary condition. These arrays have an entry
167  // corresponding to each boundary attribute. Placing a '1' in entry i
168  // marks attribute i+1 as being active, '0' is inactive.
169  Array<int> nbc_bdr(mesh->bdr_attributes.Max());
170  Array<int> rbc_bdr(mesh->bdr_attributes.Max());
171  Array<int> dbc_bdr(mesh->bdr_attributes.Max());
172 
173  nbc_bdr = 0; nbc_bdr[0] = 1;
174  rbc_bdr = 0; rbc_bdr[1] = 1;
175  dbc_bdr = 0; dbc_bdr[2] = 1;
176 
177  Array<int> ess_tdof_list(0);
178  if (h1 && mesh->bdr_attributes.Size())
179  {
180  // For a continuous basis the linear system must be modified to enforce an
181  // essential (Dirichlet) boundary condition. In the DG case this is not
182  // necessary as the boundary condition will only be enforced weakly.
183  fespace.GetEssentialTrueDofs(dbc_bdr, ess_tdof_list);
184  }
185 
186  // 5. Setup the various coefficients needed for the Laplace operator and the
187  // various boundary conditions. In general these coefficients could be
188  // functions of position but here we use only constants.
189  ConstantCoefficient matCoef(mat_val);
190  ConstantCoefficient dbcCoef(dbc_val);
191  ConstantCoefficient nbcCoef(nbc_val);
192  ConstantCoefficient rbcACoef(rbc_a_val);
193  ConstantCoefficient rbcBCoef(rbc_b_val);
194 
195  // Since the n.Grad(u) terms arise by integrating -Div(m Grad(u)) by parts we
196  // must introduce the coefficient 'm' into the boundary conditions.
197  // Therefore, in the case of the Neumann BC, we actually enforce m n.Grad(u)
198  // = m g rather than simply n.Grad(u) = g.
199  ProductCoefficient m_nbcCoef(matCoef, nbcCoef);
200  ProductCoefficient m_rbcACoef(matCoef, rbcACoef);
201  ProductCoefficient m_rbcBCoef(matCoef, rbcBCoef);
202 
203  // 6. Define the solution vector u as a finite element grid function
204  // corresponding to fespace. Initialize u with initial guess of zero.
205  GridFunction u(&fespace);
206  u = 0.0;
207 
208  // 7. Set up the bilinear form a(.,.) on the finite element space
209  // corresponding to the Laplacian operator -Delta, by adding the Diffusion
210  // domain integrator.
211  BilinearForm a(&fespace);
213  if (h1)
214  {
215  // Add a Mass integrator on the Robin boundary
216  a.AddBoundaryIntegrator(new MassIntegrator(m_rbcACoef), rbc_bdr);
217  }
218  else
219  {
220  // Add the interfacial portion of the Laplace operator
222  sigma, kappa));
223 
224  // Counteract the n.Grad(u) term on the Dirichlet portion of the boundary
225  a.AddBdrFaceIntegrator(new DGDiffusionIntegrator(matCoef, sigma, kappa),
226  dbc_bdr);
227 
228  // Augment the n.Grad(u) term with a*u on the Robin portion of boundary
230  rbc_bdr);
231  }
232  a.Assemble();
233 
234  // 8. Assemble the linear form for the right hand side vector.
235  LinearForm b(&fespace);
236 
237  if (h1)
238  {
239  // Set the Dirichlet values in the solution vector
240  u.ProjectBdrCoefficient(dbcCoef, dbc_bdr);
241 
242  // Add the desired value for n.Grad(u) on the Neumann boundary
243  b.AddBoundaryIntegrator(new BoundaryLFIntegrator(m_nbcCoef), nbc_bdr);
244 
245  // Add the desired value for n.Grad(u) + a*u on the Robin boundary
246  b.AddBoundaryIntegrator(new BoundaryLFIntegrator(m_rbcBCoef), rbc_bdr);
247  }
248  else
249  {
250  // Add the desired value for the Dirichlet boundary
251  b.AddBdrFaceIntegrator(new DGDirichletLFIntegrator(dbcCoef, matCoef,
252  sigma, kappa),
253  dbc_bdr);
254 
255  // Add the desired value for n.Grad(u) on the Neumann boundary
257  nbc_bdr);
258 
259  // Add the desired value for n.Grad(u) + a*u on the Robin boundary
260  b.AddBdrFaceIntegrator(new BoundaryLFIntegrator(m_rbcBCoef),
261  rbc_bdr);
262  }
263  b.Assemble();
264 
265  // 9. Construct the linear system.
266  OperatorPtr A;
267  Vector B, X;
268  a.FormLinearSystem(ess_tdof_list, u, b, A, X, B);
269 
270 #ifndef MFEM_USE_SUITESPARSE
271  // 10. Define a simple symmetric Gauss-Seidel preconditioner and use it to
272  // solve the system AX=B with PCG in the symmetric case, and GMRES in the
273  // non-symmetric one.
274  {
275  GSSmoother M((SparseMatrix&)(*A));
276  if (sigma == -1.0)
277  {
278  PCG(*A, M, B, X, 1, 500, 1e-12, 0.0);
279  }
280  else
281  {
282  GMRES(*A, M, B, X, 1, 500, 10, 1e-12, 0.0);
283  }
284  }
285 #else
286  // 11. If MFEM was compiled with SuiteSparse, use UMFPACK to solve the
287  // system.
288  UMFPackSolver umf_solver;
289  umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
290  umf_solver.SetOperator(*A);
291  umf_solver.Mult(B, X);
292 #endif
293 
294  // 12. Recover the grid function corresponding to U. This is the local finite
295  // element solution.
296  a.RecoverFEMSolution(X, b, u);
297 
298  // 13. Compute the various boundary integrals.
299  mfem::out << endl
300  << "Verifying boundary conditions" << endl
301  << "=============================" << endl;
302  {
303  // Integrate the solution on the Dirichlet boundary and compare to the
304  // expected value.
305  double error, avg = IntegrateBC(u, dbc_bdr, 0.0, 1.0, dbc_val, error);
306 
307  bool hom_dbc = (dbc_val == 0.0);
308  error /= hom_dbc ? 1.0 : fabs(dbc_val);
309  mfem::out << "Average of solution on Gamma_dbc:\t"
310  << avg << ", \t"
311  << (hom_dbc ? "absolute" : "relative")
312  << " error " << error << endl;
313  }
314  {
315  // Integrate n.Grad(u) on the inhomogeneous Neumann boundary and compare
316  // to the expected value.
317  double error, avg = IntegrateBC(u, nbc_bdr, 1.0, 0.0, nbc_val, error);
318 
319  bool hom_nbc = (nbc_val == 0.0);
320  error /= hom_nbc ? 1.0 : fabs(nbc_val);
321  mfem::out << "Average of n.Grad(u) on Gamma_nbc:\t"
322  << avg << ", \t"
323  << (hom_nbc ? "absolute" : "relative")
324  << " error " << error << endl;
325  }
326  {
327  // Integrate n.Grad(u) on the homogeneous Neumann boundary and compare to
328  // the expected value of zero.
329  Array<int> nbc0_bdr(mesh->bdr_attributes.Max());
330  nbc0_bdr = 0;
331  nbc0_bdr[3] = 1;
332 
333  double error, avg = IntegrateBC(u, nbc0_bdr, 1.0, 0.0, 0.0, error);
334 
335  bool hom_nbc = true;
336  mfem::out << "Average of n.Grad(u) on Gamma_nbc0:\t"
337  << avg << ", \t"
338  << (hom_nbc ? "absolute" : "relative")
339  << " error " << error << endl;
340  }
341  {
342  // Integrate n.Grad(u) + a * u on the Robin boundary and compare to the
343  // expected value.
344  double error;
345  double avg = IntegrateBC(u, rbc_bdr, 1.0, rbc_a_val, rbc_b_val, error);
346 
347  bool hom_rbc = (rbc_b_val == 0.0);
348  error /= hom_rbc ? 1.0 : fabs(rbc_b_val);
349  mfem::out << "Average of n.Grad(u)+a*u on Gamma_rbc:\t"
350  << avg << ", \t"
351  << (hom_rbc ? "absolute" : "relative")
352  << " error " << error << endl;
353  }
354 
355  // 14. Save the refined mesh and the solution. This output can be viewed
356  // later using GLVis: "glvis -m refined.mesh -g sol.gf".
357  {
358  ofstream mesh_ofs("refined.mesh");
359  mesh_ofs.precision(8);
360  mesh->Print(mesh_ofs);
361  ofstream sol_ofs("sol.gf");
362  sol_ofs.precision(8);
363  u.Save(sol_ofs);
364  }
365 
366  // 15. Send the solution by socket to a GLVis server.
367  if (visualization)
368  {
369  string title_str = h1 ? "H1" : "DG";
370  char vishost[] = "localhost";
371  int visport = 19916;
372  socketstream sol_sock(vishost, visport);
373  sol_sock.precision(8);
374  sol_sock << "solution\n" << *mesh << u
375  << "window_title '" << title_str << " Solution'"
376  << " keys 'mmc'" << flush;
377  }
378 
379  // 16. Free the used memory.
380  delete fec;
381  delete mesh;
382 
383  return 0;
384 }
385 
386 void quad_trans(double u, double v, double &x, double &y, bool log = false)
387 {
388  double a = a_; // Radius of disc
389 
390  double d = 4.0 * a * (M_SQRT2 - 2.0 * a) * (1.0 - 2.0 * v);
391 
392  double v0 = (1.0 + M_SQRT2) * (M_SQRT2 * a - 2.0 * v) *
393  ((4.0 - 3 * M_SQRT2) * a +
394  (8.0 * (M_SQRT2 - 1.0) * a - 2.0) * v) / d;
395 
396  double r = 2.0 * ((M_SQRT2 - 1.0) * a * a * (1.0 - 4.0 *v) +
397  2.0 * (1.0 + M_SQRT2 *
398  (1.0 + 2.0 * (2.0 * a - M_SQRT2 - 1.0) * a)) * v * v
399  ) / d;
400 
401  double t = asin(v / r) * u / v;
402  if (log)
403  {
404  mfem::out << "u, v, r, v0, t "
405  << u << " " << v << " " << r << " " << v0 << " " << t
406  << endl;
407  }
408  x = r * sin(t);
409  y = r * cos(t) - v0;
410 }
411 
412 void trans(const Vector &u, Vector &x)
413 {
414  double tol = 1e-4;
415 
416  if (u[1] > 0.5 - tol || u[1] < -0.5 + tol)
417  {
418  x = u;
419  return;
420  }
421  if (u[0] > 1.0 - tol || u[0] < -1.0 + tol || fabs(u[0]) < tol)
422  {
423  x = u;
424  return;
425  }
426 
427  if (u[0] > 0.0)
428  {
429  if (u[1] > fabs(u[0] - 0.5))
430  {
431  quad_trans(u[0] - 0.5, u[1], x[0], x[1]);
432  x[0] += 0.5;
433  return;
434  }
435  if (u[1] < -fabs(u[0] - 0.5))
436  {
437  quad_trans(u[0] - 0.5, -u[1], x[0], x[1]);
438  x[0] += 0.5;
439  x[1] *= -1.0;
440  return;
441  }
442  if (u[0] - 0.5 > fabs(u[1]))
443  {
444  quad_trans(u[1], u[0] - 0.5, x[1], x[0]);
445  x[0] += 0.5;
446  return;
447  }
448  if (u[0] - 0.5 < -fabs(u[1]))
449  {
450  quad_trans(u[1], 0.5 - u[0], x[1], x[0]);
451  x[0] *= -1.0;
452  x[0] += 0.5;
453  return;
454  }
455  }
456  else
457  {
458  if (u[1] > fabs(u[0] + 0.5))
459  {
460  quad_trans(u[0] + 0.5, u[1], x[0], x[1]);
461  x[0] -= 0.5;
462  return;
463  }
464  if (u[1] < -fabs(u[0] + 0.5))
465  {
466  quad_trans(u[0] + 0.5, -u[1], x[0], x[1]);
467  x[0] -= 0.5;
468  x[1] *= -1.0;
469  return;
470  }
471  if (u[0] + 0.5 > fabs(u[1]))
472  {
473  quad_trans(u[1], u[0] + 0.5, x[1], x[0]);
474  x[0] -= 0.5;
475  return;
476  }
477  if (u[0] + 0.5 < -fabs(u[1]))
478  {
479  quad_trans(u[1], -0.5 - u[0], x[1], x[0]);
480  x[0] *= -1.0;
481  x[0] -= 0.5;
482  return;
483  }
484  }
485  x = u;
486 }
487 
489 {
490  Mesh * mesh = new Mesh(2, 29, 16, 24, 2);
491 
492  int vi[4];
493 
494  for (int i=0; i<2; i++)
495  {
496  int o = 13 * i;
497  vi[0] = o + 0; vi[1] = o + 3; vi[2] = o + 4; vi[3] = o + 1;
498  mesh->AddQuad(vi);
499 
500  vi[0] = o + 1; vi[1] = o + 4; vi[2] = o + 5; vi[3] = o + 2;
501  mesh->AddQuad(vi);
502 
503  vi[0] = o + 5; vi[1] = o + 8; vi[2] = o + 9; vi[3] = o + 2;
504  mesh->AddQuad(vi);
505 
506  vi[0] = o + 8; vi[1] = o + 12; vi[2] = o + 15; vi[3] = o + 9;
507  mesh->AddQuad(vi);
508 
509  vi[0] = o + 11; vi[1] = o + 14; vi[2] = o + 15; vi[3] = o + 12;
510  mesh->AddQuad(vi);
511 
512  vi[0] = o + 10; vi[1] = o + 13; vi[2] = o + 14; vi[3] = o + 11;
513  mesh->AddQuad(vi);
514 
515  vi[0] = o + 6; vi[1] = o + 13; vi[2] = o + 10; vi[3] = o + 7;
516  mesh->AddQuad(vi);
517 
518  vi[0] = o + 0; vi[1] = o + 6; vi[2] = o + 7; vi[3] = o + 3;
519  mesh->AddQuad(vi);
520  }
521 
522  vi[0] = 0; vi[1] = 6; mesh->AddBdrSegment(vi, 1);
523  vi[0] = 6; vi[1] = 13; mesh->AddBdrSegment(vi, 1);
524  vi[0] = 13; vi[1] = 19; mesh->AddBdrSegment(vi, 1);
525  vi[0] = 19; vi[1] = 26; mesh->AddBdrSegment(vi, 1);
526 
527  vi[0] = 28; vi[1] = 22; mesh->AddBdrSegment(vi, 2);
528  vi[0] = 22; vi[1] = 15; mesh->AddBdrSegment(vi, 2);
529  vi[0] = 15; vi[1] = 9; mesh->AddBdrSegment(vi, 2);
530  vi[0] = 9; vi[1] = 2; mesh->AddBdrSegment(vi, 2);
531 
532  for (int i=0; i<2; i++)
533  {
534  int o = 13 * i;
535  vi[0] = o + 7; vi[1] = o + 3; mesh->AddBdrSegment(vi, 3 + i);
536  vi[0] = o + 10; vi[1] = o + 7; mesh->AddBdrSegment(vi, 3 + i);
537  vi[0] = o + 11; vi[1] = o + 10; mesh->AddBdrSegment(vi, 3 + i);
538  vi[0] = o + 12; vi[1] = o + 11; mesh->AddBdrSegment(vi, 3 + i);
539  vi[0] = o + 8; vi[1] = o + 12; mesh->AddBdrSegment(vi, 3 + i);
540  vi[0] = o + 5; vi[1] = o + 8; mesh->AddBdrSegment(vi, 3 + i);
541  vi[0] = o + 4; vi[1] = o + 5; mesh->AddBdrSegment(vi, 3 + i);
542  vi[0] = o + 3; vi[1] = o + 4; mesh->AddBdrSegment(vi, 3 + i);
543  }
544 
545  double d[2];
546  double a = a_ / M_SQRT2;
547 
548  d[0] = -1.0; d[1] = -0.5; mesh->AddVertex(d);
549  d[0] = -1.0; d[1] = 0.0; mesh->AddVertex(d);
550  d[0] = -1.0; d[1] = 0.5; mesh->AddVertex(d);
551 
552  d[0] = -0.5 - a; d[1] = -a; mesh->AddVertex(d);
553  d[0] = -0.5 - a; d[1] = 0.0; mesh->AddVertex(d);
554  d[0] = -0.5 - a; d[1] = a; mesh->AddVertex(d);
555 
556  d[0] = -0.5; d[1] = -0.5; mesh->AddVertex(d);
557  d[0] = -0.5; d[1] = -a; mesh->AddVertex(d);
558  d[0] = -0.5; d[1] = a; mesh->AddVertex(d);
559  d[0] = -0.5; d[1] = 0.5; mesh->AddVertex(d);
560 
561  d[0] = -0.5 + a; d[1] = -a; mesh->AddVertex(d);
562  d[0] = -0.5 + a; d[1] = 0.0; mesh->AddVertex(d);
563  d[0] = -0.5 + a; d[1] = a; mesh->AddVertex(d);
564 
565  d[0] = 0.0; d[1] = -0.5; mesh->AddVertex(d);
566  d[0] = 0.0; d[1] = 0.0; mesh->AddVertex(d);
567  d[0] = 0.0; d[1] = 0.5; mesh->AddVertex(d);
568 
569  d[0] = 0.5 - a; d[1] = -a; mesh->AddVertex(d);
570  d[0] = 0.5 - a; d[1] = 0.0; mesh->AddVertex(d);
571  d[0] = 0.5 - a; d[1] = a; mesh->AddVertex(d);
572 
573  d[0] = 0.5; d[1] = -0.5; mesh->AddVertex(d);
574  d[0] = 0.5; d[1] = -a; mesh->AddVertex(d);
575  d[0] = 0.5; d[1] = a; mesh->AddVertex(d);
576  d[0] = 0.5; 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] = 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  mesh->FinalizeTopology();
587 
588  mesh->SetCurvature(1, true);
589 
590  // Stitch the ends of the stack together
591  {
592  Array<int> v2v(mesh->GetNV());
593  for (int i = 0; i < v2v.Size() - 3; i++)
594  {
595  v2v[i] = i;
596  }
597  // identify vertices on the narrow ends of the rectangle
598  v2v[v2v.Size() - 3] = 0;
599  v2v[v2v.Size() - 2] = 1;
600  v2v[v2v.Size() - 1] = 2;
601 
602  // renumber elements
603  for (int i = 0; i < mesh->GetNE(); i++)
604  {
605  Element *el = mesh->GetElement(i);
606  int *v = el->GetVertices();
607  int nv = el->GetNVertices();
608  for (int j = 0; j < nv; j++)
609  {
610  v[j] = v2v[v[j]];
611  }
612  }
613  // renumber boundary elements
614  for (int i = 0; i < mesh->GetNBE(); i++)
615  {
616  Element *el = mesh->GetBdrElement(i);
617  int *v = el->GetVertices();
618  int nv = el->GetNVertices();
619  for (int j = 0; j < nv; j++)
620  {
621  v[j] = v2v[v[j]];
622  }
623  }
624  mesh->RemoveUnusedVertices();
625  mesh->RemoveInternalBoundaries();
626  }
627  mesh->SetCurvature(3, true);
628 
629  for (int l = 0; l < ref; l++)
630  {
631  mesh->UniformRefinement();
632  }
633 
634  mesh->Transform(trans);
635 
636  return mesh;
637 }
638 
639 double IntegrateBC(const GridFunction &x, const Array<int> &bdr,
640  double alpha, double beta, double gamma,
641  double &error)
642 {
643  double nrm = 0.0;
644  double avg = 0.0;
645  error = 0.0;
646 
647  const bool a_is_zero = alpha == 0.0;
648  const bool b_is_zero = beta == 0.0;
649 
650  const FiniteElementSpace &fes = *x.FESpace();
651  MFEM_ASSERT(fes.GetVDim() == 1, "");
652  Mesh &mesh = *fes.GetMesh();
653  Vector shape, loc_dofs, w_nor;
654  DenseMatrix dshape;
655  Array<int> dof_ids;
656  for (int i = 0; i < mesh.GetNBE(); i++)
657  {
658  if (bdr[mesh.GetBdrAttribute(i)-1] == 0) { continue; }
659 
660  FaceElementTransformations *FTr = mesh.GetBdrFaceTransformations(i);
661  if (FTr == nullptr) { continue; }
662 
663  const FiniteElement &fe = *fes.GetFE(FTr->Elem1No);
664  MFEM_ASSERT(fe.GetMapType() == FiniteElement::VALUE, "");
665  const int int_order = 2*fe.GetOrder() + 3;
666  const IntegrationRule &ir = IntRules.Get(FTr->FaceGeom, int_order);
667 
668  fes.GetElementDofs(FTr->Elem1No, dof_ids);
669  x.GetSubVector(dof_ids, loc_dofs);
670  if (!a_is_zero)
671  {
672  const int sdim = FTr->Face->GetSpaceDim();
673  w_nor.SetSize(sdim);
674  dshape.SetSize(fe.GetDof(), sdim);
675  }
676  if (!b_is_zero)
677  {
678  shape.SetSize(fe.GetDof());
679  }
680  for (int j = 0; j < ir.GetNPoints(); j++)
681  {
682  const IntegrationPoint &ip = ir.IntPoint(j);
683  IntegrationPoint eip;
684  FTr->Loc1.Transform(ip, eip);
685  FTr->Face->SetIntPoint(&ip);
686  double face_weight = FTr->Face->Weight();
687  double val = 0.0;
688  if (!a_is_zero)
689  {
690  FTr->Elem1->SetIntPoint(&eip);
691  fe.CalcPhysDShape(*FTr->Elem1, dshape);
692  CalcOrtho(FTr->Face->Jacobian(), w_nor);
693  val += alpha * dshape.InnerProduct(w_nor, loc_dofs) / face_weight;
694  }
695  if (!b_is_zero)
696  {
697  fe.CalcShape(eip, shape);
698  val += beta * (shape * loc_dofs);
699  }
700 
701  // Measure the length of the boundary
702  nrm += ip.weight * face_weight;
703 
704  // Integrate alpha * n.Grad(x) + beta * x
705  avg += val * ip.weight * face_weight;
706 
707  // Integrate |alpha * n.Grad(x) + beta * x - gamma|^2
708  val -= gamma;
709  error += (val*val) * ip.weight * face_weight;
710  }
711  }
712 
713  // Normalize by the length of the boundary
714  if (std::abs(nrm) > 0.0)
715  {
716  error /= nrm;
717  avg /= nrm;
718  }
719 
720  // Compute l2 norm of the error in the boundary condition (negative
721  // quadrature weights may produce negative 'error')
722  error = (error >= 0.0) ? sqrt(error) : -sqrt(-error);
723 
724  // Return the average value of alpha * n.Grad(x) + beta * x
725  return avg;
726 }
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_base.hpp:235
FDualNumber< tbase > asin(const FDualNumber< tbase > &f)
asin([dual number])
Definition: fdual.hpp:439
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
void trans(const Vector &u, Vector &x)
Definition: ex27.cpp:412
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
int AddQuad(int v1, int v2, int v3, int v4, int attr=1)
Definition: mesh.cpp:1662
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:923
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
A coefficient that is constant across space and time.
Definition: coefficient.hpp:78
virtual void GetVertices(Array< int > &v) const =0
Returns element&#39;s vertices.
void n4Vec(const Vector &x, Vector &n)
Definition: ex27.cpp:69
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:521
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:931
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:792
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:102
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:284
void quad_trans(double u, double v, double &x, double &y, bool log=false)
Definition: ex27.cpp:386
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
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition: vector.cpp:534
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 RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
Recover the solution of a linear system formed with FormLinearSystem().
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
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:11551
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Get a list of essential true dofs, ess_tdof_list, corresponding to the boundary attributes marked in ...
Definition: fespace.cpp:564
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:928
double kappa
Definition: ex24.cpp:54
int GetMapType() const
Returns the FiniteElement::MapType of the element describing how reference functions are mapped to ph...
Definition: fe_base.hpp:349
void CalcOrtho(const DenseMatrix &J, Vector &n)
Definition: densemat.cpp:2288
void RemoveInternalBoundaries()
Definition: mesh.cpp:11707
Data type for Gauss-Seidel smoother of sparse matrix.
Direct sparse solver using UMFPACK.
Definition: solvers.hpp:1025
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:3619
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:1601
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
Data type sparse matrix.
Definition: sparsemat.hpp:46
constexpr char vishost[]
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:433
IntegrationPointTransformation Loc1
Definition: eltrans.hpp:524
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9737
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:119
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:5312
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
int AddBdrSegment(int v1, int v2, int attr=1)
Definition: mesh.cpp:1811
double Weight()
Return the weight of the Jacobian matrix of the transformation at the currently set IntegrationPoint...
Definition: eltrans.hpp:131
const Element * GetElement(int i) const
Definition: mesh.hpp:1033
void PCG(const Operator &A, Solver &B, const Vector &b, Vector &x, int print_iter, int max_num_iter, double RTOLERANCE, double ATOLERANCE)
Preconditioned conjugate gradient method. (tolerances are squared)
Definition: solvers.cpp:904
void AddBdrFaceIntegrator(LinearFormIntegrator *lfi)
Adds new Boundary Face Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:83
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:566
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:191
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:557
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:652
void AddBoundaryIntegrator(LinearFormIntegrator *lfi)
Adds new Boundary Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:70
int Dimension() const
Definition: mesh.hpp:999
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
virtual DofTransformation * GetElementDofs(int elem, Array< int > &dofs) const
Returns indices of degrees of freedom of element &#39;elem&#39;.
Definition: fespace.cpp:2680
FDualNumber< tbase > cos(const FDualNumber< tbase > &f)
cos([dual number])
Definition: fdual.hpp:471
L2_FECollection DG_FECollection
Declare an alternative name for L2_FECollection = DG_FECollection.
Definition: fe_coll.hpp:334
double IntegrateBC(const GridFunction &sol, const Array< int > &bdr_marker, double alpha, double beta, double gamma, double &error)
Definition: ex27.cpp:639
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:270
double Control[UMFPACK_CONTROL]
Definition: solvers.hpp:1036
FDualNumber< tbase > sin(const FDualNumber< tbase > &f)
sin([dual number])
Definition: fdual.hpp:572
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:88
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:323
void RemoveUnusedVertices()
Remove unused vertices and rebuild mesh connectivity.
Definition: mesh.cpp:11600
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:2878
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 Transform(const IntegrationPoint &, IntegrationPoint &)
Definition: eltrans.cpp:546
FDualNumber< tbase > log(const FDualNumber< tbase > &f)
log([dual number])
Definition: fdual.hpp:515
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:925
A &quot;square matrix&quot; operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
Class for integration point with weight.
Definition: intrules.hpp:25
virtual void Print(std::ostream &os=mfem::out) const
Definition: mesh.hpp:1646
FDualNumber< tbase > sqrt(const FDualNumber< tbase > &f)
sqrt([dual number])
Definition: fdual.hpp:600
void AddBoundaryIntegrator(BilinearFormIntegrator *bfi)
Adds new Boundary Integrator. Assumes ownership of bfi.
int dim
Definition: ex24.cpp:53
void AddInteriorFaceIntegrator(BilinearFormIntegrator *bfi)
Adds new interior Face Integrator. Assumes ownership of bfi.
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
Mesh * GenerateSerialMesh(int ref)
Definition: ex27.cpp:488
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th element in t...
Definition: fespace.cpp:2783
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
ElementTransformation * Elem1
Definition: eltrans.hpp:522
ElementTransformation * Face
Definition: eltrans.hpp:523
virtual int GetNVertices() const =0
RefCoord t[3]
const double alpha
Definition: ex15.cpp:369
Vector data type.
Definition: vector.hpp:60
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:216
Vector with associated FE space and LinearFormIntegrators.
Definition: linearform.hpp:23
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
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
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:379
Abstract data type element.
Definition: element.hpp:28
int main()
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:3179
const Element * GetBdrElement(int i) const
Definition: mesh.hpp:1037
int GMRES(const Operator &A, Vector &x, const Vector &b, Solver &M, int &max_iter, int m, double &tol, double atol, int printit)
GMRES method. (tolerances are squared)
Definition: solvers.cpp:1335
void AddBdrFaceIntegrator(BilinearFormIntegrator *bfi)
Adds new boundary Face Integrator. Assumes ownership of bfi.
void ProjectBdrCoefficient(Coefficient &coeff, Array< int > &attr)
Project a Coefficient on the GridFunction, modifying only DOFs on the boundary associated with the bo...
Definition: gridfunc.hpp:447
virtual void SetOperator(const Operator &op)
Factorize the given Operator op which must be a SparseMatrix.
Definition: solvers.cpp:3084
double sigma(const Vector &x)
Definition: maxwell.cpp:102
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150