MFEM  v3.3 Finite element discretization library
ex14.cpp
Go to the documentation of this file.
1 // MFEM Example 14
2 //
3 // Compile with: make ex14
4 //
5 // Sample runs: ex14 -m ../data/inline-quad.mesh -o 0
6 // ex14 -m ../data/star.mesh -r 4 -o 2
7 // ex14 -m ../data/escher.mesh -s 1
8 // ex14 -m ../data/fichera.mesh -s 1 -k 1
9 // ex14 -m ../data/square-disc-p2.vtk -r 3 -o 2
10 // ex14 -m ../data/square-disc-p3.mesh -r 2 -o 3
11 // ex14 -m ../data/square-disc-nurbs.mesh -o 1
12 // ex14 -m ../data/disc-nurbs.mesh -r 3 -o 2 -s 1 -k 0
13 // ex14 -m ../data/pipe-nurbs.mesh -o 1
14 // ex14 -m ../data/inline-segment.mesh -r 5
15 // ex14 -m ../data/amr-quad.mesh -r 3
16 // ex14 -m ../data/amr-hex.mesh
17 // ex14 -m ../data/fichera-amr.mesh
18 //
19 // Description: This example code demonstrates the use of MFEM to define a
20 // discontinuous Galerkin (DG) finite element discretization of
21 // the Laplace problem -Delta u = 1 with homogeneous Dirichlet
22 // boundary conditions. Finite element spaces of any order,
23 // including zero on regular grids, are supported. The example
24 // highlights the use of discontinuous spaces and DG-specific face
25 // integrators.
26 //
27 // We recommend viewing examples 1 and 9 before viewing this
28 // example.
29
30 #include "mfem.hpp"
31 #include <fstream>
32 #include <iostream>
33
34 using namespace std;
35 using namespace mfem;
36
37 int main(int argc, char *argv[])
38 {
39  // 1. Parse command-line options.
40  const char *mesh_file = "../data/star.mesh";
41  int ref_levels = -1;
42  int order = 1;
43  double sigma = -1.0;
44  double kappa = -1.0;
45  bool visualization = 1;
46
47  OptionsParser args(argc, argv);
48  args.AddOption(&mesh_file, "-m", "--mesh",
49  "Mesh file to use.");
50  args.AddOption(&ref_levels, "-r", "--refine",
51  "Number of times to refine the mesh uniformly, -1 for auto.");
52  args.AddOption(&order, "-o", "--order",
53  "Finite element order (polynomial degree) >= 0.");
54  args.AddOption(&sigma, "-s", "--sigma",
55  "One of the two DG penalty parameters, typically +1/-1."
56  " See the documentation of class DGDiffusionIntegrator.");
57  args.AddOption(&kappa, "-k", "--kappa",
58  "One of the two DG penalty parameters, should be positive."
59  " Negative values are replaced with (order+1)^2.");
60  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
61  "--no-visualization",
62  "Enable or disable GLVis visualization.");
63  args.Parse();
64  if (!args.Good())
65  {
66  args.PrintUsage(cout);
67  return 1;
68  }
69  if (kappa < 0)
70  {
71  kappa = (order+1)*(order+1);
72  }
73  args.PrintOptions(cout);
74
75  // 2. Read the mesh from the given mesh file. We can handle triangular,
76  // quadrilateral, tetrahedral and hexahedral meshes with the same code.
77  // NURBS meshes are projected to second order meshes.
78  Mesh *mesh = new Mesh(mesh_file, 1, 1);
79  int dim = mesh->Dimension();
80
81  // 3. Refine the mesh to increase the resolution. In this example we do
82  // 'ref_levels' of uniform refinement. By default, or if ref_levels < 0,
83  // we choose it to be the largest number that gives a final mesh with no
84  // more than 50,000 elements.
85  {
86  if (ref_levels < 0)
87  {
88  ref_levels = (int)floor(log(50000./mesh->GetNE())/log(2.)/dim);
89  }
90  for (int l = 0; l < ref_levels; l++)
91  {
92  mesh->UniformRefinement();
93  }
94  }
95  if (mesh->NURBSext)
96  {
97  mesh->SetCurvature(max(order, 1));
98  }
99
100  // 4. Define a finite element space on the mesh. Here we use discontinuous
101  // finite elements of the specified order >= 0.
102  FiniteElementCollection *fec = new DG_FECollection(order, dim);
103  FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec);
104  cout << "Number of unknowns: " << fespace->GetVSize() << endl;
105
106  // 5. Set up the linear form b(.) which corresponds to the right-hand side of
107  // the FEM linear system.
108  LinearForm *b = new LinearForm(fespace);
109  ConstantCoefficient one(1.0);
110  ConstantCoefficient zero(0.0);
113  new DGDirichletLFIntegrator(zero, one, sigma, kappa));
114  b->Assemble();
115
116  // 6. Define the solution vector x as a finite element grid function
117  // corresponding to fespace. Initialize x with initial guess of zero.
118  GridFunction x(fespace);
119  x = 0.0;
120
121  // 7. Set up the bilinear form a(.,.) on the finite element space
122  // corresponding to the Laplacian operator -Delta, by adding the Diffusion
123  // domain integrator and the interior and boundary DG face integrators.
124  // Note that boundary conditions are imposed weakly in the form, so there
125  // is no need for dof elimination. After assembly and finalizing we
126  // extract the corresponding sparse matrix A.
127  BilinearForm *a = new BilinearForm(fespace);
129  a->AddInteriorFaceIntegrator(new DGDiffusionIntegrator(one, sigma, kappa));
130  a->AddBdrFaceIntegrator(new DGDiffusionIntegrator(one, sigma, kappa));
131  a->Assemble();
132  a->Finalize();
133  const SparseMatrix &A = a->SpMat();
134
135 #ifndef MFEM_USE_SUITESPARSE
136  // 8. Define a simple symmetric Gauss-Seidel preconditioner and use it to
137  // solve the system Ax=b with PCG in the symmetric case, and GMRES in the
138  // non-symmetric one.
139  GSSmoother M(A);
140  if (sigma == -1.0)
141  {
142  PCG(A, M, *b, x, 1, 500, 1e-12, 0.0);
143  }
144  else
145  {
146  GMRES(A, M, *b, x, 1, 500, 10, 1e-12, 0.0);
147  }
148 #else
149  // 8. If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
150  UMFPackSolver umf_solver;
151  umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
152  umf_solver.SetOperator(A);
153  umf_solver.Mult(*b, x);
154 #endif
155
156  // 9. Save the refined mesh and the solution. This output can be viewed later
157  // using GLVis: "glvis -m refined.mesh -g sol.gf".
158  ofstream mesh_ofs("refined.mesh");
159  mesh_ofs.precision(8);
160  mesh->Print(mesh_ofs);
161  ofstream sol_ofs("sol.gf");
162  sol_ofs.precision(8);
163  x.Save(sol_ofs);
164
165  // 10. Send the solution by socket to a GLVis server.
166  if (visualization)
167  {
168  char vishost[] = "localhost";
169  int visport = 19916;
170  socketstream sol_sock(vishost, visport);
171  sol_sock.precision(8);
172  sol_sock << "solution\n" << *mesh << x << flush;
173  }
174
175  // 11. Free the used memory.
176  delete a;
177  delete b;
178  delete fespace;
179  delete fec;
180  delete mesh;
181
182  return 0;
183 }
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:47
int GetVSize() const
Definition: fespace.hpp:163
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:27
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
Subclass constant coefficient.
Definition: coefficient.hpp:57
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:42
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:584
int main(int argc, char *argv[])
Data type for Gauss-Seidel smoother of sparse matrix.
Direct sparse solver using UMFPACK.
Definition: solvers.hpp:337
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:2185
int dim
Definition: ex3.cpp:47
Data type sparse matrix.
Definition: sparsemat.hpp:38
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:6684
void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:3278
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:456
void AddBdrFaceIntegrator(LinearFormIntegrator *lfi)
Adds new Boundary Face Integrator.
Definition: linearform.cpp:29
int Dimension() const
Definition: mesh.hpp:611
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
L2_FECollection DG_FECollection
Definition: fe_coll.hpp:236
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator.
Definition: linearform.cpp:19
double Control[UMFPACK_CONTROL]
Definition: solvers.hpp:348
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)
Definition: optparser.hpp:74
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:144
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
void AddInteriorFaceIntegrator(BilinearFormIntegrator *bfi)
Adds new interior Face Integrator.
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator.
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
double kappa
Definition: ex3.cpp:46
Class for linear form - Vector with associated FE space and LFIntegrators.
Definition: linearform.hpp:23
const SparseMatrix & SpMat() const
Returns a reference to the sparse matrix.
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:1774
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:843
void AddBdrFaceIntegrator(BilinearFormIntegrator *bfi)
Adds new boundary Face Integrator.
virtual void SetOperator(const Operator &op)
Factorize the given Operator op which must be a SparseMatrix.
Definition: solvers.cpp:1674
virtual void Print(std::ostream &out=std::cout) const
Definition: mesh.hpp:988
bool Good() const
Definition: optparser.hpp:120