MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
ex27.cpp
Go to the documentation of this file.
1// MFEM Example 27
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// | | \ |
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 GridFunction &sol, const Array<int> &bdr_marker,
78 real_t &error);
79
80int main(int argc, char *argv[])
81{
82 // 1. Parse command-line options.
83 int ser_ref_levels = 2;
84 int order = 1;
85 real_t sigma = -1.0;
86 real_t kappa = -1.0;
87 bool h1 = true;
88 bool visualization = true;
89
90 real_t mat_val = 1.0;
91 real_t dbc_val = 0.0;
92 real_t nbc_val = 1.0;
93 real_t rbc_a_val = 1.0; // du/dn + a * u = b
94 real_t 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.");
100 "Finite element order (polynomial degree) or -1 for"
101 " isoparametric space.");
103 "One of the two DG penalty parameters, typically +1/-1."
104 " See the documentation of class DGDiffusionIntegrator.");
106 "One of the two DG penalty parameters, should be positive."
107 " Negative values are replaced with (order+1)^2.");
109 "Number of times to refine the mesh uniformly in serial.");
111 "Constant value for material coefficient "
112 "in the Laplace operator.");
114 "Constant value for Dirichlet Boundary Condition.");
116 "Constant value for Neumann Boundary Condition.");
118 "Constant 'a' value for Robin Boundary Condition: "
119 "du/dn + a * u = b.");
121 "Constant 'b' value for Robin Boundary Condition: "
122 "du/dn + a * u = b.");
124 "Radius of holes in the mesh.");
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 }
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) :
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
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
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
244
245 // Add the desired value for n.Grad(u) + a*u on the Robin boundary
247 }
248 else
249 {
250 // Add the desired value for the Dirichlet boundary
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
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 real_t 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 real_t 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 real_t 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 real_t error;
345 real_t 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
386void quad_trans(real_t u, real_t v, real_t &x, real_t &y, bool log = false)
387{
388 real_t a = a_; // Radius of disc
389
390 real_t d = 4.0 * a * (M_SQRT2 - 2.0 * a) * (1.0 - 2.0 * v);
391
392 real_t 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 real_t 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 real_t 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
412void trans(const Vector &u, Vector &x)
413{
414 real_t 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;
499
500 vi[0] = o + 1; vi[1] = o + 4; vi[2] = o + 5; vi[3] = o + 2;
502
503 vi[0] = o + 5; vi[1] = o + 8; vi[2] = o + 9; vi[3] = o + 2;
505
506 vi[0] = o + 8; vi[1] = o + 12; vi[2] = o + 15; vi[3] = o + 9;
508
509 vi[0] = o + 11; vi[1] = o + 14; vi[2] = o + 15; vi[3] = o + 12;
511
512 vi[0] = o + 10; vi[1] = o + 13; vi[2] = o + 14; vi[3] = o + 11;
514
515 vi[0] = o + 6; vi[1] = o + 13; vi[2] = o + 10; vi[3] = o + 7;
517
518 vi[0] = o + 0; vi[1] = o + 6; vi[2] = o + 7; vi[3] = o + 3;
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 real_t d[2];
546 real_t 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();
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
640 real_t alpha, real_t beta, real_t gamma,
641 real_t &error)
642{
643 real_t nrm = 0.0;
644 real_t 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
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);
684 FTr->Loc1.Transform(ip, eip);
685 FTr->Face->SetIntPoint(&ip);
686 real_t face_weight = FTr->Face->Weight();
687 real_t 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}
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
A "square matrix" operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
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
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:220
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition fespace.hpp:716
DofTransformation * GetElementDofs(int elem, Array< int > &dofs) const
Returns indices of degrees of freedom of element 'elem'. The returned indices are offsets into an ldo...
Definition fespace.cpp:2846
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1) const
Get a list of essential true dofs, ess_tdof_list, corresponding to the boundary attributes marked in ...
Definition fespace.cpp:588
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th element in t...
Definition fespace.cpp:3168
Mesh * GetMesh() const
Returns the mesh.
Definition fespace.hpp:559
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
Data type for Gauss-Seidel smoother of sparse matrix.
Class for grid function - Vector with associated FE space.
Definition gridfunc.hpp:31
FiniteElementSpace * FESpace()
Definition gridfunc.hpp:696
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
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.
Vector with associated FE space and LinearFormIntegrators.
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
virtual void Print(std::ostream &os=mfem::out, const std::string &comments="") const
Definition mesh.hpp:2288
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
Pointer to an Operator of a specified type.
Definition handle.hpp:34
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.
Scalar coefficient defined as the product of two scalar coefficients or a scalar and a scalar coeffic...
Data type sparse matrix.
Definition sparsemat.hpp:51
Direct sparse solver using UMFPACK.
Definition solvers.hpp:1096
virtual void SetOperator(const Operator &op)
Factorize the given Operator op which must be a SparseMatrix.
Definition solvers.cpp:3102
real_t Control[UMFPACK_CONTROL]
Definition solvers.hpp:1106
virtual void Mult(const Vector &b, Vector &x) const
Direct solution of the linear system using UMFPACK.
Definition solvers.cpp:3197
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 ex27.cpp:69
void trans(const Vector &u, Vector &x)
Definition ex27.cpp:412
real_t IntegrateBC(const GridFunction &sol, const Array< int > &bdr_marker, real_t alpha, real_t beta, real_t gamma, real_t &error)
Definition ex27.cpp:639
void quad_trans(real_t u, real_t v, real_t &x, real_t &y, bool log=false)
Definition ex27.cpp:386
Mesh * GenerateSerialMesh(int ref)
Definition ex27.cpp:488
int main()
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
int GMRES(const Operator &A, Vector &x, const Vector &b, Solver &M, int &max_iter, int m, real_t &tol, real_t atol, int printit)
GMRES method. (tolerances are squared)
Definition solvers.cpp:1342
void PCG(const Operator &A, Solver &B, const Vector &b, Vector &x, int print_iter, int max_num_iter, real_t RTOLERANCE, real_t ATOLERANCE)
Preconditioned conjugate gradient method. (tolerances are squared)
Definition solvers.cpp:913
L2_FECollection DG_FECollection
Declare an alternative name for L2_FECollection = DG_FECollection.
Definition fe_coll.hpp:382
float real_t
Definition config.hpp:43