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