MFEM  v3.0
ex3.cpp
Go to the documentation of this file.
1 // MFEM Example 3
2 //
3 // Compile with: make ex3
4 //
5 // Sample runs: ex3 -m ../data/beam-tet.mesh
6 // ex3 -m ../data/beam-hex.mesh
7 // ex3 -m ../data/escher.mesh
8 // ex3 -m ../data/fichera.mesh
9 // ex3 -m ../data/fichera-q2.vtk
10 // ex3 -m ../data/fichera-q3.mesh
11 // ex3 -m ../data/beam-hex-nurbs.mesh
12 //
13 // Description: This example code solves a simple 3D electromagnetic diffusion
14 // problem corresponding to the second order definite Maxwell
15 // equation curl curl E + E = f with boundary condition
16 // E x n = <given tangential field>. Here, we use a given exact
17 // solution E and compute the corresponding r.h.s. f.
18 // We discretize with Nedelec finite elements.
19 //
20 // The example demonstrates the use of H(curl) finite element
21 // spaces with the curl-curl and the (vector finite element) mass
22 // bilinear form, as well as the computation of discretization
23 // error when the exact solution is known.
24 //
25 // We recommend viewing examples 1-2 before viewing this example.
26
27 #include "mfem.hpp"
28 #include <fstream>
29 #include <iostream>
30
31 using namespace std;
32 using namespace mfem;
33
34 // Exact solution, E, and r.h.s., f. See below for implementation.
35 void E_exact(const Vector &, Vector &);
36 void f_exact(const Vector &, Vector &);
37
38 int main(int argc, char *argv[])
39 {
40  // 1. Parse command-line options.
41  const char *mesh_file = "../data/beam-tet.mesh";
42  int order = 1;
43  bool visualization = 1;
44
45  OptionsParser args(argc, argv);
47  "Mesh file to use.");
49  "Finite element order (polynomial degree).");
51  "--no-visualization",
52  "Enable or disable GLVis visualization.");
53  args.Parse();
54  if (!args.Good())
55  {
56  args.PrintUsage(cout);
57  return 1;
58  }
59  args.PrintOptions(cout);
60
61  // 2. Read the mesh from the given mesh file. We can handle triangular,
62  // quadrilateral, tetrahedral, hexahedral, surface and volume meshes with
63  // the same code.
64  Mesh *mesh;
65  ifstream imesh(mesh_file);
66  if (!imesh)
67  {
68  cerr << "\nCan not open mesh file: " << mesh_file << '\n' << endl;
69  return 2;
70  }
71  mesh = new Mesh(imesh, 1, 1);
72  imesh.close();
73  int dim = mesh->Dimension();
74  if (dim != 3)
75  {
76  cerr << "\nThis example requires a 3D mesh\n" << endl;
77  return 3;
78  }
79
80  // 3. Refine the mesh to increase the resolution. In this example we do
81  // 'ref_levels' of uniform refinement. We choose 'ref_levels' to be the
82  // largest number that gives a final mesh with no more than 50,000
83  // elements.
84  {
85  int ref_levels =
86  (int)floor(log(50000./mesh->GetNE())/log(2.)/dim);
87  for (int l = 0; l < ref_levels; l++)
88  mesh->UniformRefinement();
89  }
90  mesh->ReorientTetMesh();
91
92  // 4. Define a finite element space on the mesh. Here we use the lowest order
93  // Nedelec finite elements, but we can easily switch to higher-order
94  // spaces by changing the value of p.
95  FiniteElementCollection *fec = new ND_FECollection(order, dim);
96  FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec);
97  cout << "Number of unknowns: " << fespace->GetVSize() << endl;
98
99  // 5. Set up the linear form b(.) which corresponds to the right-hand side
100  // of the FEM linear system, which in this case is (f,phi_i) where f is
101  // given by the function f_exact and phi_i are the basis functions in the
102  // finite element fespace.
104  LinearForm *b = new LinearForm(fespace);
106  b->Assemble();
107
108  // 6. Define the solution vector x as a finite element grid function
109  // corresponding to fespace. Initialize x by projecting the exact
110  // solution. Note that only values from the boundary edges will be used
111  // when eliminating the non-homogeneous boundary condition to modify the
112  // r.h.s. vector b.
113  GridFunction x(fespace);
115  x.ProjectCoefficient(E);
116
117  // 7. Set up the bilinear form corresponding to the EM diffusion operator
118  // curl muinv curl + sigma I, by adding the curl-curl and the mass domain
119  // integrators and finally imposing the non-homogeneous Dirichlet boundary
120  // conditions. The boundary conditions are implemented by marking all the
121  // boundary attributes from the mesh as essential (Dirichlet). After
122  // assembly and finalizing we extract the corresponding sparse matrix A.
123  Coefficient *muinv = new ConstantCoefficient(1.0);
124  Coefficient *sigma = new ConstantCoefficient(1.0);
125  BilinearForm *a = new BilinearForm(fespace);
128  a->Assemble();
129  Array<int> ess_bdr(mesh->bdr_attributes.Max());
130  ess_bdr = 1;
131  a->EliminateEssentialBC(ess_bdr, x, *b);
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.
138  GSSmoother M(A);
139  x = 0.0;
140  PCG(A, M, *b, x, 1, 500, 1e-12, 0.0);
141 #else
142  // 8. If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
143  UMFPackSolver umf_solver;
144  umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
145  umf_solver.SetOperator(A);
146  umf_solver.Mult(*b, x);
147 #endif
148
149  // 9. Compute and print the L^2 norm of the error.
150  cout << "\n|| E_h - E ||_{L^2} = " << x.ComputeL2Error(E) << '\n' << endl;
151
152  // 10. Save the refined mesh and the solution. This output can be viewed
153  // later using GLVis: "glvis -m refined.mesh -g sol.gf".
154  {
155  ofstream mesh_ofs("refined.mesh");
156  mesh_ofs.precision(8);
157  mesh->Print(mesh_ofs);
158  ofstream sol_ofs("sol.gf");
159  sol_ofs.precision(8);
160  x.Save(sol_ofs);
161  }
162
163  // 11. Send the solution by socket to a GLVis server.
164  if (visualization)
165  {
166  char vishost[] = "localhost";
167  int visport = 19916;
168  socketstream sol_sock(vishost, visport);
169  sol_sock.precision(8);
170  sol_sock << "solution\n" << *mesh << x << flush;
171  }
172
173  // 12. Free the used memory.
174  delete a;
175  delete sigma;
176  delete muinv;
177  delete b;
178  delete fespace;
179  delete fec;
180  delete mesh;
181
182  return 0;
183 }
184
185 // A parameter for the exact solution.
186 const double kappa = M_PI;
187
188 void E_exact(const Vector &x, Vector &E)
189 {
190  E(0) = sin(kappa * x(1));
191  E(1) = sin(kappa * x(2));
192  E(2) = sin(kappa * x(0));
193 }
194
195 void f_exact(const Vector &x, Vector &f)
196 {
197  f(0) = (1. + kappa * kappa) * sin(kappa * x(1));
198  f(1) = (1. + kappa * kappa) * sin(kappa * x(2));
199  f(2) = (1. + kappa * kappa) * sin(kappa * x(0));
200 }
int GetVSize() const
Definition: fespace.hpp:151
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:26
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
Subclass constant coefficient.
Definition: coefficient.hpp:57
Integrator for (curl u, curl v) for Nedelec elements.
Definition: bilininteg.hpp:388
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:396
Data type for Gauss-Seidel smoother of sparse matrix.
Direct sparse solver using UMFPACK.
Definition: solvers.hpp:329
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:1797
double ComputeL2Error(Coefficient &exsol, const IntegrationRule *irs[]=NULL) const
Definition: gridfunc.hpp:166
Data type sparse matrix.
Definition: sparsemat.hpp:38
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:6225
T Max() const
Definition: array.cpp:78
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:424
int Dimension() const
Definition: mesh.hpp:417
virtual void ReorientTetMesh()
Definition: mesh.cpp:4118
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:385
const double kappa
Definition: ex3.cpp:186
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:7136
Definition: linearform.cpp:19
Array< int > bdr_attributes
Definition: mesh.hpp:305
int main(int argc, char *argv[])
Definition: ex1.cpp:39
double Control[UMFPACK_CONTROL]
Definition: solvers.hpp:340
Abstract finite element space.
Definition: fespace.hpp:61
Base class Coefficient that may optionally depend on time.
Definition: coefficient.hpp:31
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
void f_exact(const Vector &, Vector &)
Definition: ex3.cpp:195
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:266
void ProjectCoefficient(Coefficient &coeff)
Definition: gridfunc.cpp:967
void E_exact(const Vector &, Vector &)
Definition: ex3.cpp:188
for VectorFiniteElements (Nedelec, Raviart-Thomas)
Definition: lininteg.hpp:171
Arbitrary order H(curl)-conforming Nedelec finite elements.
Definition: fe_coll.hpp:158
Vector data type.
Definition: vector.hpp:29
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 martix.
void EliminateEssentialBC(Array< int > &bdr_attr_is_ess, Vector &sol, Vector &rhs, int d=0)
Integrator for (Q u, v) for VectorFiniteElements.
Definition: bilininteg.hpp:434
virtual void Mult(const Vector &b, Vector &x) const
Operator application.
Definition: solvers.cpp:1597
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: solvers.cpp:1505
bool Good() const
Definition: optparser.hpp:120