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