MFEM  v3.0
ex3p.cpp
Go to the documentation of this file.
1 // MFEM Example 3 - Parallel Version
2 //
3 // Compile with: make ex3p
4 //
5 // Sample runs: mpirun -np 4 ex3p -m ../data/beam-tet.mesh
6 // mpirun -np 4 ex3p -m ../data/beam-hex.mesh
7 // mpirun -np 4 ex3p -m ../data/escher.mesh
8 // mpirun -np 4 ex3p -m ../data/fichera.mesh
9 // mpirun -np 4 ex3p -m ../data/fichera-q2.vtk
10 // mpirun -np 4 ex3p -m ../data/fichera-q3.mesh
11 // mpirun -np 4 ex3p -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. Initialize MPI.
41  int num_procs, myid;
42  MPI_Init(&argc, &argv);
43  MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
44  MPI_Comm_rank(MPI_COMM_WORLD, &myid);
45
46  // 2. Parse command-line options.
47  const char *mesh_file = "../data/beam-tet.mesh";
48  int order = 1;
49  bool visualization = 1;
50
51  OptionsParser args(argc, argv);
53  "Mesh file to use.");
55  "Finite element order (polynomial degree).");
57  "--no-visualization",
58  "Enable or disable GLVis visualization.");
59  args.Parse();
60  if (!args.Good())
61  {
62  if (myid == 0)
63  args.PrintUsage(cout);
64  MPI_Finalize();
65  return 1;
66  }
67  if (myid == 0)
68  args.PrintOptions(cout);
69
70  // 3. Read the (serial) mesh from the given mesh file on all processors. We
71  // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
72  // and volume meshes with the same code.
73  Mesh *mesh;
74  ifstream imesh(mesh_file);
75  if (!imesh)
76  {
77  if (myid == 0)
78  cerr << "\nCan not open mesh file: " << mesh_file << '\n' << endl;
79  MPI_Finalize();
80  return 2;
81  }
82  mesh = new Mesh(imesh, 1, 1);
83  imesh.close();
84  int dim = mesh->Dimension();
85  if (dim != 3)
86  {
87  if (myid == 0)
88  cerr << "\nThis example requires a 3D mesh\n" << endl;
89  MPI_Finalize();
90  return 3;
91  }
92
93  // 4. Refine the serial mesh on all processors to increase the resolution. In
94  // this example we do 'ref_levels' of uniform refinement. We choose
95  // 'ref_levels' to be the largest number that gives a final mesh with no
96  // more than 1,000 elements.
97  {
98  int ref_levels =
99  (int)floor(log(1000./mesh->GetNE())/log(2.)/dim);
100  for (int l = 0; l < ref_levels; l++)
101  mesh->UniformRefinement();
102  }
103
104  // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
105  // this mesh further in parallel to increase the resolution. Once the
106  // parallel mesh is defined, the serial mesh can be deleted. Tetrahedral
107  // meshes need to be reoriented before we can define high-order Nedelec
108  // spaces on them.
109  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
110  delete mesh;
111  {
112  int par_ref_levels = 2;
113  for (int l = 0; l < par_ref_levels; l++)
114  pmesh->UniformRefinement();
115  }
116  pmesh->ReorientTetMesh();
117
118  // 6. Define a parallel finite element space on the parallel mesh. Here we
119  // use the lowest order Nedelec finite elements, but we can easily switch
120  // to higher-order spaces by changing the value of p.
121  FiniteElementCollection *fec = new ND_FECollection(order, dim);
122  ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
123  int size = fespace->GlobalTrueVSize();
124  if (myid == 0)
125  cout << "Number of unknowns: " << size << endl;
126
127  // 7. Set up the parallel linear form b(.) which corresponds to the
128  // right-hand side of the FEM linear system, which in this case is
129  // (f,phi_i) where f is given by the function f_exact and phi_i are the
130  // basis functions in the finite element fespace.
132  ParLinearForm *b = new ParLinearForm(fespace);
134  b->Assemble();
135
136  // 8. Define the solution vector x as a parallel finite element grid function
137  // corresponding to fespace. Initialize x by projecting the exact
138  // solution. Note that only values from the boundary edges will be used
139  // when eliminating the non-homogeneous boundary condition to modify the
140  // r.h.s. vector b.
141  ParGridFunction x(fespace);
143  x.ProjectCoefficient(E);
144
145  // 9. Set up the parallel bilinear form corresponding to the EM diffusion
146  // operator curl muinv curl + sigma I, by adding the curl-curl and the
147  // mass domain integrators and finally imposing non-homogeneous Dirichlet
148  // boundary conditions. The boundary conditions are implemented by
149  // marking all the boundary attributes from the mesh as essential
150  // (Dirichlet). After serial and parallel assembly we extract the
151  // parallel matrix A.
152  Coefficient *muinv = new ConstantCoefficient(1.0);
153  Coefficient *sigma = new ConstantCoefficient(1.0);
154  ParBilinearForm *a = new ParBilinearForm(fespace);
157  a->Assemble();
158  Array<int> ess_bdr(pmesh->bdr_attributes.Max());
159  ess_bdr = 1;
160  a->EliminateEssentialBC(ess_bdr, x, *b);
161  a->Finalize();
162
163  // 10. Define the parallel (hypre) matrix and vectors representing a(.,.),
164  // b(.) and the finite element approximation.
168  *X = 0.0;
169
170  delete a;
171  delete sigma;
172  delete muinv;
173  delete b;
174
175  // 11. Define and apply a parallel PCG solver for AX=B with the AMS
176  // preconditioner from hypre.
177  HypreSolver *ams = new HypreAMS(*A, fespace);
178  HyprePCG *pcg = new HyprePCG(*A);
179  pcg->SetTol(1e-12);
180  pcg->SetMaxIter(500);
181  pcg->SetPrintLevel(2);
182  pcg->SetPreconditioner(*ams);
183  pcg->Mult(*B, *X);
184
185  // 12. Extract the parallel grid function corresponding to the finite element
186  // approximation X. This is the local solution on each processor.
187  x = *X;
188
189  // 13. Compute and print the L^2 norm of the error.
190  {
191  double err = x.ComputeL2Error(E);
192  if (myid == 0)
193  cout << "\n|| E_h - E ||_{L^2} = " << err << '\n' << endl;
194  }
195
196  // 14. Save the refined mesh and the solution in parallel. This output can
197  // be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
198  {
199  ostringstream mesh_name, sol_name;
200  mesh_name << "mesh." << setfill('0') << setw(6) << myid;
201  sol_name << "sol." << setfill('0') << setw(6) << myid;
202
203  ofstream mesh_ofs(mesh_name.str().c_str());
204  mesh_ofs.precision(8);
205  pmesh->Print(mesh_ofs);
206
207  ofstream sol_ofs(sol_name.str().c_str());
208  sol_ofs.precision(8);
209  x.Save(sol_ofs);
210  }
211
212  // 15. Send the solution by socket to a GLVis server.
213  if (visualization)
214  {
215  char vishost[] = "localhost";
216  int visport = 19916;
217  socketstream sol_sock(vishost, visport);
218  sol_sock << "parallel " << num_procs << " " << myid << "\n";
219  sol_sock.precision(8);
220  sol_sock << "solution\n" << *pmesh << x << flush;
221  }
222
223  // 16. Free the used memory.
224  delete pcg;
225  delete ams;
226  delete X;
227  delete B;
228  delete A;
229  delete fespace;
230  delete fec;
231  delete pmesh;
232
233  MPI_Finalize();
234
235  return 0;
236 }
237
238 // A parameter for the exact solution.
239 const double kappa = M_PI;
240
241 void E_exact(const Vector &x, Vector &E)
242 {
243  E(0) = sin(kappa * x(1));
244  E(1) = sin(kappa * x(2));
245  E(2) = sin(kappa * x(0));
246 }
247
248 void f_exact(const Vector &x, Vector &f)
249 {
250  f(0) = (1. + kappa * kappa) * sin(kappa * x(1));
251  f(1) = (1. + kappa * kappa) * sin(kappa * x(2));
252  f(2) = (1. + kappa * kappa) * sin(kappa * x(0));
253 }
void SetTol(double tol)
Definition: hypre.cpp:1339
The Auxiliary-space Maxwell Solver in hypre.
Definition: hypre.hpp:554
double ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL) const
Definition: pgridfunc.hpp:128
Subclass constant coefficient.
Definition: coefficient.hpp:57
virtual void ReorientTetMesh()
See the remarks for the serial version in mesh.hpp.
Definition: pmesh.cpp:1200
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
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:250
Abstract parallel finite element space.
Definition: pfespace.hpp:28
void ProjectCoefficient(Coefficient &coeff)
Definition: pgridfunc.cpp:229
void SetPrintLevel(int print_lvl)
Definition: hypre.cpp:1354
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
Class for parallel linear form.
Definition: plinearform.hpp:26
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:6225
T Max() const
Definition: array.cpp:78
void Assemble(int skip_zeros=1)
Assemble the local matrix.
int Dimension() const
Definition: mesh.hpp:417
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:385
void SetMaxIter(int max_iter)
Definition: hypre.cpp:1344
const double kappa
Definition: ex3.cpp:186
Wrapper for hypre&#39;s parallel vector class.
Definition: hypre.hpp:38
Definition: linearform.cpp:19
Array< int > bdr_attributes
Definition: mesh.hpp:305
int main(int argc, char *argv[])
Definition: ex1.cpp:39
PCG solver in hypre.
Definition: hypre.hpp:381
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 SetPreconditioner(HypreSolver &precond)
Set the hypre solver to be used as a preconditioner.
Definition: hypre.cpp:1360
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:266
void E_exact(const Vector &, Vector &)
Definition: ex3.cpp:188
Class for parallel bilinear form.
Abstract class for hypre&#39;s solvers and preconditioners.
Definition: hypre.hpp:345
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 parallel grid function.
Definition: pgridfunc.hpp:31
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:103
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Solve Ax=b with hypre&#39;s PCG.
Definition: hypre.cpp:1377
Class for parallel meshes.
Definition: pmesh.hpp:27
void ParallelAverage(Vector &tv) const
Returns the vector averaged on the true dofs.
Definition: pgridfunc.cpp:110
void EliminateEssentialBC(Array< int > &bdr_attr_is_ess, Vector &sol, Vector &rhs, int d=0)
void ParallelAssemble(Vector &tv)
Assemble the vector on the true dofs, i.e. P^t v.
Definition: plinearform.cpp:34
Integrator for (Q u, v) for VectorFiniteElements.
Definition: bilininteg.hpp:434
virtual void Print(std::ostream &out=std::cout) const
Definition: pmesh.cpp:2444
bool Good() const
Definition: optparser.hpp:120