MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
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// | | \ |
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
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
63using namespace std;
64using namespace mfem;
65
66static real_t a_ = 0.2;
67
68// Normal to hole with boundary attribute 4
69void n4Vec(const Vector &x, Vector &n) { n = x; n[0] -= 0.5; n /= -n.Norml2(); }
70
71Mesh * 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.
76real_t IntegrateBC(const ParGridFunction &sol, const Array<int> &bdr_marker,
78 real_t &error);
79
80int main(int argc, char *argv[])
81{
82 // 1. Initialize MPI and HYPRE.
83 Mpi::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 real_t sigma = -1.0;
92 real_t kappa = -1.0;
93 bool h1 = true;
94 bool visualization = true;
95
96 real_t mat_val = 1.0;
97 real_t dbc_val = 0.0;
98 real_t nbc_val = 1.0;
99 real_t rbc_a_val = 1.0; // du/dn + a * u = b
100 real_t 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.");
106 "Finite element order (polynomial degree) or -1 for"
107 " isoparametric space.");
109 "One of the two DG penalty parameters, typically +1/-1."
110 " See the documentation of class DGDiffusionIntegrator.");
112 "One of the two DG penalty parameters, should be positive."
113 " Negative values are replaced with (order+1)^2.");
115 "Number of times to refine the mesh uniformly in serial.");
117 "Number of times to refine the mesh uniformly in parallel.");
119 "Constant value for material coefficient "
120 "in the Laplace operator.");
122 "Constant value for Dirichlet Boundary Condition.");
124 "Constant value for Neumann Boundary Condition.");
126 "Constant 'a' value for Robin Boundary Condition: "
127 "du/dn + a * u = b.");
129 "Constant 'b' value for Robin Boundary Condition: "
130 "du/dn + a * u = b.");
132 "Radius of holes in the mesh.");
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 }
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);
231 if (h1)
232 {
233 // Add a Mass integrator on the Robin boundary
235 }
236 else
237 {
238 // Add the interfacial portion of the Laplace operator
240 sigma, kappa));
241
242 // Counteract the n.Grad(u) term on the Dirichlet portion of the boundary
244 dbc_bdr);
245
246 // Augment the n.Grad(u) term with a*u on the Robin portion of boundary
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
262
263 // Add the desired value for n.Grad(u) + a*u on the Robin boundary
265 }
266 else
267 {
268 // Add the desired value for the Dirichlet boundary
270 sigma, kappa),
271 dbc_bdr);
272
273 // Add the desired value for n.Grad(u) on the Neumann boundary
275 nbc_bdr);
276
277 // Add the desired value for n.Grad(u) + a*u on the Robin boundary
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 real_t 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 real_t 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 real_t 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 real_t 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
412void quad_trans(real_t u, real_t v, real_t &x, real_t &y, bool log = false)
413{
414 real_t a = a_; // Radius of disc
415
416 real_t d = 4.0 * a * (M_SQRT2 - 2.0 * a) * (1.0 - 2.0 * v);
417
418 real_t 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 real_t 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 real_t 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
438void trans(const Vector &u, Vector &x)
439{
440 real_t 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;
525
526 vi[0] = o + 1; vi[1] = o + 4; vi[2] = o + 5; vi[3] = o + 2;
528
529 vi[0] = o + 5; vi[1] = o + 8; vi[2] = o + 9; vi[3] = o + 2;
531
532 vi[0] = o + 8; vi[1] = o + 12; vi[2] = o + 15; vi[3] = o + 9;
534
535 vi[0] = o + 11; vi[1] = o + 14; vi[2] = o + 15; vi[3] = o + 12;
537
538 vi[0] = o + 10; vi[1] = o + 13; vi[2] = o + 14; vi[3] = o + 11;
540
541 vi[0] = o + 6; vi[1] = o + 13; vi[2] = o + 10; vi[3] = o + 7;
543
544 vi[0] = o + 0; vi[1] = o + 6; vi[2] = o + 7; vi[3] = o + 3;
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 real_t d[2];
572 real_t 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();
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
666 real_t alpha, real_t beta, real_t gamma,
667 real_t &glb_err)
668{
669 real_t loc_vals[3];
670 real_t &nrm = loc_vals[0];
671 real_t &avg = loc_vals[1];
672 real_t &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
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);
715 FTr->Loc1.Transform(ip, eip);
716 FTr->Face->SetIntPoint(&ip);
717 real_t face_weight = FTr->Face->Weight();
718 real_t 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 real_t glb_vals[3];
745 MPI_Allreduce(loc_vals, glb_vals, 3, MPITypeMap<real_t>::mpi_type,
746 MPI_SUM, fes.GetComm());
747
748 real_t glb_nrm = glb_vals[0];
749 real_t glb_avg = glb_vals[1];
750 glb_err = glb_vals[2];
751
752 // Normalize by the length of the boundary
753 if (std::abs(glb_nrm) > 0.0)
754 {
755 glb_err /= glb_nrm;
756 glb_avg /= glb_nrm;
757 }
758
759 // Compute l2 norm of the error in the boundary condition (negative
760 // quadrature weights may produce negative 'error')
761 glb_err = (glb_err >= 0.0) ? sqrt(glb_err) : -sqrt(-glb_err);
762
763 // Return the average value of alpha * n.Grad(x) + beta * x
764 return glb_avg;
765}
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
int Size() const
Return the logical size of the array.
Definition array.hpp:144
Class for boundary integration .
Definition lininteg.hpp:180
A coefficient that is constant across space and time.
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
real_t InnerProduct(const real_t *x, const real_t *y) const
Compute y^t A x.
Definition densemat.cpp:383
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition densemat.hpp:105
virtual int GetSpaceDim() const =0
Get the dimension of the target (physical) space.
real_t Weight()
Return the weight of the Jacobian matrix of the transformation at the currently set IntegrationPoint....
Definition eltrans.hpp:131
const DenseMatrix & Jacobian()
Return the Jacobian matrix of the transformation at the currently set IntegrationPoint,...
Definition eltrans.hpp:119
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition eltrans.hpp:93
Abstract data type element.
Definition element.hpp:29
virtual void GetVertices(Array< int > &v) const =0
Get the indices defining the vertices.
virtual int GetNVertices() const =0
A specialized ElementTransformation class representing a face and its two neighboring elements.
Definition eltrans.hpp:484
ElementTransformation * Elem1
Definition eltrans.hpp:525
ElementTransformation * Face
Definition eltrans.hpp:526
IntegrationPointTransformation Loc1
Definition eltrans.hpp:527
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition fe_coll.hpp:27
int GetVDim() const
Returns vector dimension.
Definition fespace.hpp:706
Abstract class for all finite elements.
Definition fe_base.hpp:239
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition fe_base.hpp:333
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
int GetMapType() const
Returns the FiniteElement::MapType of the element describing how reference functions are mapped to ph...
Definition fe_base.hpp:355
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...
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition fe_base.hpp:329
GMRES method.
Definition solvers.hpp:547
void SetKDim(int dim)
Set the number of iteration to perform between restarts, default is 50.
Definition solvers.hpp:559
virtual void Mult(const Vector &b, Vector &x) const
Iterative solution of the linear system using the GMRES method.
Definition solvers.cpp:980
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
The BoomerAMG solver in hypre.
Definition hypre.hpp:1691
PCG solver in hypre.
Definition hypre.hpp:1275
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition hypre.cpp:4114
void SetPrintLevel(int print_lvl)
Definition hypre.cpp:4156
void SetPreconditioner(HypreSolver &precond)
Set the hypre solver to be used as a preconditioner.
Definition hypre.cpp:4161
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Solve Ax=b with hypre's PCG.
Definition hypre.cpp:4184
void SetMaxIter(int max_iter)
Definition hypre.cpp:4146
void SetTol(real_t tol)
Definition hypre.cpp:4136
Abstract class for hypre's solvers and preconditioners.
Definition hypre.hpp:1162
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
Definition hypre.hpp:74
void Transform(const IntegrationPoint &, IntegrationPoint &)
Definition eltrans.cpp:543
Class for integration point with weight.
Definition intrules.hpp:35
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:100
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition intrules.hpp:256
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition intrules.hpp:259
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
virtual void SetOperator(const Operator &op) override
Also calls SetOperator for the preconditioner if there is one.
Definition solvers.cpp:179
void SetRelTol(real_t rtol)
Definition solvers.hpp:209
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition solvers.cpp:173
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition solvers.cpp:71
void SetMaxIter(int max_it)
Definition solvers.hpp:211
void SetAbsTol(real_t atol)
Definition solvers.hpp:210
Mesh data type.
Definition mesh.hpp:56
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition mesh.hpp:282
Adds a quadrilateral to the mesh given by 4 vertices v1 through v4.
Definition mesh.cpp:1743
int GetBdrAttribute(int i) const
Return the attribute of boundary element i.
Definition mesh.hpp:1339
const Element * GetElement(int i) const
Return pointer to the i'th element object.
Definition mesh.hpp:1283
void FinalizeTopology(bool generate_bdr=true)
Finalize the construction of the secondary topology (connectivity) data of a Mesh.
Definition mesh.cpp:3135
int AddVertex(real_t x, real_t y=0.0, real_t z=0.0)
Definition mesh.cpp:1658
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1226
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
const Element * GetBdrElement(int i) const
Return pointer to the i'th boundary element object.
Definition mesh.hpp:1298
int AddBdrSegment(int v1, int v2, int attr=1)
Definition mesh.cpp:2035
FaceElementTransformations * GetBdrFaceTransformations(int BdrElemNo)
Builds the transformation defining the given boundary face.
Definition mesh.cpp:1099
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
Definition mesh.hpp:1223
int GetNBE() const
Returns number of boundary elements.
Definition mesh.hpp:1229
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Set the curvature of the mesh nodes using the given polynomial degree.
Definition mesh.cpp:6211
void Transform(void(*f)(const Vector &, Vector &))
Definition mesh.cpp:12821
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
void RemoveInternalBoundaries()
Definition mesh.cpp:12979
void RemoveUnusedVertices()
Remove unused vertices and rebuild mesh connectivity.
Definition mesh.cpp:12872
static bool Root()
Return true if the rank in MPI_COMM_WORLD is zero.
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
static int WorldSize()
Return the size of MPI_COMM_WORLD.
static void Init(int &argc, char **&argv, int required=default_thread_required, int *provided=nullptr)
Singleton creation with Mpi::Init(argc, argv).
Pointer to an Operator of a specified type.
Definition handle.hpp:34
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
Reconstruct a solution vector x (e.g. a GridFunction) from the solution X of a constrained linear sys...
Definition operator.cpp:148
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void PrintUsage(std::ostream &out) const
Print the usage message.
void PrintOptions(std::ostream &out) const
Print the options.
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
Add a boolean option and set 'var' to receive the value. Enable/disable tags are used to set the bool...
Definition optparser.hpp:82
bool Good() const
Return true if the command line options were parsed successfully.
void Disable()
Disable output.
Definition globals.hpp:54
Class for parallel bilinear form.
Abstract parallel finite element space.
Definition pfespace.hpp:29
MPI_Comm GetComm() const
Definition pfespace.hpp:273
void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1) const override
HYPRE_BigInt GlobalTrueVSize() const
Definition pfespace.hpp:285
void GetElementDofs(int i, Array< int > &dofs, DofTransformation &doftrans) const override
The same as GetElementDofs(), but with a user-allocated DofTransformation object. doftrans must be al...
Definition pfespace.cpp:468
ParMesh * GetParMesh() const
Definition pfespace.hpp:277
const FiniteElement * GetFE(int i) const override
Definition pfespace.cpp:534
Class for parallel grid function.
Definition pgridfunc.hpp:33
ParFiniteElementSpace * ParFESpace() const
Class for parallel linear form.
Class for parallel meshes.
Definition pmesh.hpp:34
void Print(std::ostream &out=mfem::out, const std::string &comments="") const override
Definition pmesh.cpp:4801
Scalar coefficient defined as the product of two scalar coefficients or a scalar and a scalar coeffic...
Vector data type.
Definition vector.hpp:80
real_t Norml2() const
Returns the l2 norm of the vector.
Definition vector.cpp:832
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:538
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition vector.cpp:578
Vector beta
real_t sigma(const Vector &x)
Definition maxwell.cpp:102
const real_t alpha
Definition ex15.cpp:369
real_t kappa
Definition ex24.cpp:54
int dim
Definition ex24.cpp:53
void n4Vec(const Vector &x, Vector &n)
Definition ex27p.cpp:69
real_t IntegrateBC(const ParGridFunction &sol, const Array< int > &bdr_marker, real_t alpha, real_t beta, real_t gamma, real_t &error)
Definition ex27p.cpp:665
void trans(const Vector &u, Vector &x)
Definition ex27p.cpp:438
void quad_trans(real_t u, real_t v, real_t &x, real_t &y, bool log=false)
Definition ex27p.cpp:412
Mesh * GenerateSerialMesh(int ref)
Definition ex27p.cpp:514
int main()
HYPRE_Int HYPRE_BigInt
real_t b
Definition lissajous.cpp:42
real_t a
Definition lissajous.cpp:41
const int visport
void CalcOrtho(const DenseMatrix &J, Vector &n)
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition globals.hpp:66
L2_FECollection DG_FECollection
Declare an alternative name for L2_FECollection = DG_FECollection.
Definition fe_coll.hpp:382
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
float real_t
Definition config.hpp:43