MFEM  v3.1 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);
49  "Mesh file to use.");
51  "Number of times to refine the mesh uniformly, -1 for auto.");
53  "Finite element order (polynomial degree) >= 0.");
55  "One of the two DG penalty parameters, typically +1/-1."
56  " See the documentation of class DiffusionIntegrator.");
58  "One of the two DG penalty parameters, should be positive."
59  " Negative values are replaced with (order+1)^2.");
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;
79  ifstream imesh(mesh_file);
80  if (!imesh)
81  {
82  cerr << "\nCan not open mesh file: " << mesh_file << '\n' << endl;
83  return 2;
84  }
85  mesh = new Mesh(imesh, 1, 1);
86  imesh.close();
87  int dim = mesh->Dimension();
88  if (mesh->NURBSext)
89  {
90  mesh->SetCurvature(2);
91  }
92
93  // 3. Refine the mesh to increase the resolution. In this example we do
94  // 'ref_levels' of uniform refinement. By default, or if ref_levels < 0,
95  // we choose it to be the largest number that gives a final mesh with no
96  // more than 50,000 elements.
97  {
98  if (ref_levels < 0)
99  {
100  ref_levels = (int)floor(log(50000./mesh->GetNE())/log(2.)/dim);
101  }
102  for (int l = 0; l < ref_levels; l++)
103  {
104  mesh->UniformRefinement();
105  }
106  }
107
108  // 4. Define a finite element space on the mesh. Here we use discontinuous
109  // finite elements of the specified order >= 0.
110  FiniteElementCollection *fec = new DG_FECollection(order, dim);
111  FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec);
112  cout << "Number of unknowns: " << fespace->GetVSize() << endl;
113
114  // 5. Set up the linear form b(.) which corresponds to the right-hand side of
115  // the FEM linear system.
116  LinearForm *b = new LinearForm(fespace);
117  ConstantCoefficient one(1.0);
118  ConstantCoefficient zero(0.0);
121  new DGDirichletLFIntegrator(zero, one, sigma, kappa));
122  b->Assemble();
123
124  // 6. Define the solution vector x as a finite element grid function
125  // corresponding to fespace. Initialize x with initial guess of zero.
126  GridFunction x(fespace);
127  x = 0.0;
128
129  // 7. Set up the bilinear form a(.,.) on the finite element space
130  // corresponding to the Laplacian operator -Delta, by adding the Diffusion
131  // domain integrator and the interior and boundary DG face integrators.
132  // Note that boundary conditions are imposed weakly in the form, so there
133  // is no need for dof elimination. After assembly and finalizing we
134  // extract the corresponding sparse matrix A.
135  BilinearForm *a = new BilinearForm(fespace);
139  a->Assemble();
140  a->Finalize();
141  const SparseMatrix &A = a->SpMat();
142
143 #ifndef MFEM_USE_SUITESPARSE
144  // 8. Define a simple symmetric Gauss-Seidel preconditioner and use it to
145  // solve the system Ax=b with PCG in the symmetric case, and GMRES in the
146  // non-symmetric one.
147  GSSmoother M(A);
148  if (sigma == -1.0)
149  {
150  PCG(A, M, *b, x, 1, 500, 1e-12, 0.0);
151  }
152  else
153  {
154  GMRES(A, M, *b, x, 1, 500, 10, 1e-12, 0.0);
155  }
156 #else
157  // 8. If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
158  UMFPackSolver umf_solver;
159  umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
160  umf_solver.SetOperator(A);
161  umf_solver.Mult(*b, x);
162 #endif
163
164  // 9. Save the refined mesh and the solution. This output can be viewed later
165  // using GLVis: "glvis -m refined.mesh -g sol.gf".
166  ofstream mesh_ofs("refined.mesh");
167  mesh_ofs.precision(8);
168  mesh->Print(mesh_ofs);
169  ofstream sol_ofs("sol.gf");
170  sol_ofs.precision(8);
171  x.Save(sol_ofs);
172
173  // 10. Send the solution by socket to a GLVis server.
174  if (visualization)
175  {
176  char vishost[] = "localhost";
177  int visport = 19916;
178  socketstream sol_sock(vishost, visport);
179  sol_sock.precision(8);
180  sol_sock << "solution\n" << *mesh << x << flush;
181  }
182
183  // 11. Free the used memory.
184  delete a;
185  delete b;
186  delete fespace;
187  delete fec;
188  delete mesh;
189
190  return 0;
191 }
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:47
int GetVSize() const
Definition: fespace.hpp:164
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:34
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:454
Data type for Gauss-Seidel smoother of sparse matrix.
Direct sparse solver using UMFPACK.
Definition: solvers.hpp:332
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:2138
int dim
Definition: ex3.cpp:48
Data type sparse matrix.
Definition: sparsemat.hpp:38
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:7316
void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:3736
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:440
Definition: linearform.cpp:29
int Dimension() const
Definition: mesh.hpp:475
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
L2_FECollection DG_FECollection
Definition: fe_coll.hpp:130
virtual void Print(std::ostream &out=std::cout) const
Print the mesh to the given stream using the default MFEM mesh format.
Definition: mesh.cpp:8332
Definition: linearform.cpp:19
int main(int argc, char *argv[])
Definition: ex1.cpp:45
double Control[UMFPACK_CONTROL]
Definition: solvers.hpp:343
Abstract finite element space.
Definition: fespace.hpp:62
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
Definition: mesh.hpp:142
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
double kappa
Definition: ex3.cpp:47
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.
Definition: solvers.cpp:1725
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:808
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: solvers.cpp:1625
bool Good() const
Definition: optparser.hpp:120