MFEM  v4.6.0
Finite element discretization library
ex26.cpp
Go to the documentation of this file.
1 // MFEM Example 26
2 //
3 // Compile with: make ex26
4 //
5 // Sample runs: ex26 -m ../data/star.mesh
6 // ex26 -m ../data/fichera.mesh
7 // ex26 -m ../data/beam-hex.mesh
8 //
9 // Device sample runs:
10 // ex26 -d cuda
11 // ex26 -d raja-cuda
12 // ex26 -d occa-cuda
13 // ex26 -d raja-omp
14 // ex26 -d occa-omp
15 // ex26 -d ceed-cpu
16 // ex26 -d ceed-cuda
17 // ex26 -m ../data/beam-hex.mesh -d cuda
18 //
19 // Description: This example code demonstrates the use of MFEM to define a
20 // simple finite element discretization of the Laplace problem
21 // -Delta u = 1 with homogeneous Dirichlet boundary conditions
22 // as in Example 1.
23 //
24 // It highlights on the creation of a hierarchy of discretization
25 // spaces with partial assembly and the construction of an
26 // efficient multigrid preconditioner for the iterative solver.
27 //
28 // We recommend viewing Example 1 before viewing this example.
29 
30 #include "mfem.hpp"
31 #include <fstream>
32 #include <iostream>
33 
34 using namespace std;
35 using namespace mfem;
36 
37 // Class for constructing a multigrid preconditioner for the diffusion operator.
38 // This example multigrid preconditioner class demonstrates the creation of the
39 // diffusion bilinear forms and operators using partial assembly for all spaces
40 // in the FiniteElementSpaceHierarchy. The preconditioner uses a CG solver on
41 // the coarsest level and second order Chebyshev accelerated smoothers on the
42 // other levels.
43 class DiffusionMultigrid : public GeometricMultigrid
44 {
45 private:
47 
48 public:
49  // Constructs a diffusion multigrid for the given FiniteElementSpaceHierarchy
50  // and the array of essential boundaries
51  DiffusionMultigrid(FiniteElementSpaceHierarchy& fespaces, Array<int>& ess_bdr)
52  : GeometricMultigrid(fespaces), one(1.0)
53  {
54  ConstructCoarseOperatorAndSolver(fespaces.GetFESpaceAtLevel(0), ess_bdr);
55 
56  for (int level = 1; level < fespaces.GetNumLevels(); ++level)
57  {
58  ConstructOperatorAndSmoother(fespaces.GetFESpaceAtLevel(level), ess_bdr);
59  }
60  }
61 
62 private:
63  void ConstructBilinearForm(FiniteElementSpace& fespace, Array<int>& ess_bdr)
64  {
65  BilinearForm* form = new BilinearForm(&fespace);
66  form->SetAssemblyLevel(AssemblyLevel::PARTIAL);
68  form->Assemble();
69  bfs.Append(form);
70 
71  essentialTrueDofs.Append(new Array<int>());
72  fespace.GetEssentialTrueDofs(ess_bdr, *essentialTrueDofs.Last());
73  }
74 
75  void ConstructCoarseOperatorAndSolver(FiniteElementSpace& coarse_fespace,
76  Array<int>& ess_bdr)
77  {
78  ConstructBilinearForm(coarse_fespace, ess_bdr);
79 
80  OperatorPtr opr;
81  opr.SetType(Operator::ANY_TYPE);
82  bfs.Last()->FormSystemMatrix(*essentialTrueDofs.Last(), opr);
83  opr.SetOperatorOwner(false);
84 
85  CGSolver* pcg = new CGSolver();
86  pcg->SetPrintLevel(-1);
87  pcg->SetMaxIter(200);
88  pcg->SetRelTol(sqrt(1e-4));
89  pcg->SetAbsTol(0.0);
90  pcg->SetOperator(*opr.Ptr());
91 
92  AddLevel(opr.Ptr(), pcg, true, true);
93  }
94 
95  void ConstructOperatorAndSmoother(FiniteElementSpace& fespace,
96  Array<int>& ess_bdr)
97  {
98  ConstructBilinearForm(fespace, ess_bdr);
99 
100  OperatorPtr opr;
101  opr.SetType(Operator::ANY_TYPE);
102  bfs.Last()->FormSystemMatrix(*essentialTrueDofs.Last(), opr);
103  opr.SetOperatorOwner(false);
104 
105  Vector diag(fespace.GetTrueVSize());
106  bfs.Last()->AssembleDiagonal(diag);
107 
108  Solver* smoother = new OperatorChebyshevSmoother(*opr, diag,
109  *essentialTrueDofs.Last(), 2);
110  AddLevel(opr.Ptr(), smoother, true, true);
111  }
112 };
113 
114 
115 int main(int argc, char *argv[])
116 {
117  // 1. Parse command-line options.
118  const char *mesh_file = "../data/star.mesh";
119  int geometric_refinements = 0;
120  int order_refinements = 2;
121  const char *device_config = "cpu";
122  bool visualization = true;
123 
124  OptionsParser args(argc, argv);
125  args.AddOption(&mesh_file, "-m", "--mesh",
126  "Mesh file to use.");
127  args.AddOption(&geometric_refinements, "-gr", "--geometric-refinements",
128  "Number of geometric refinements done prior to order refinements.");
129  args.AddOption(&order_refinements, "-or", "--order-refinements",
130  "Number of order refinements. Finest level in the hierarchy has order 2^{or}.");
131  args.AddOption(&device_config, "-d", "--device",
132  "Device configuration string, see Device::Configure().");
133  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
134  "--no-visualization",
135  "Enable or disable GLVis visualization.");
136  args.Parse();
137  if (!args.Good())
138  {
139  args.PrintUsage(cout);
140  return 1;
141  }
142  args.PrintOptions(cout);
143 
144  // 2. Enable hardware devices such as GPUs, and programming models such as
145  // CUDA, OCCA, RAJA and OpenMP based on command line options.
146  Device device(device_config);
147  device.Print();
148 
149  // 3. Read the mesh from the given mesh file. We can handle triangular,
150  // quadrilateral, tetrahedral, hexahedral, surface and volume meshes with
151  // the same code.
152  Mesh *mesh = new Mesh(mesh_file, 1, 1);
153  int dim = mesh->Dimension();
154 
155  // 4. Refine the mesh to increase the resolution. In this example we do
156  // 'ref_levels' of uniform refinement. We choose 'ref_levels' to be the
157  // largest number that gives a final mesh with no more than 50,000
158  // elements.
159  {
160  int ref_levels =
161  (int)floor(log(5000./mesh->GetNE())/log(2.)/dim);
162  for (int l = 0; l < ref_levels; l++)
163  {
164  mesh->UniformRefinement();
165  }
166  }
167 
168  // 5. Define a finite element space hierarchy on the mesh. Here we use
169  // continuous Lagrange finite elements. We start with order 1 on the
170  // coarse level and geometrically refine the spaces by the specified
171  // amount. Afterwards, we increase the order of the finite elements
172  // by a factor of 2 for each additional level.
174  FiniteElementSpace *coarse_fespace = new FiniteElementSpace(mesh, fec);
175  FiniteElementSpaceHierarchy fespaces(mesh, coarse_fespace, true, true);
176 
178  collections.Append(fec);
179  for (int level = 0; level < geometric_refinements; ++level)
180  {
181  fespaces.AddUniformlyRefinedLevel();
182  }
183  for (int level = 0; level < order_refinements; ++level)
184  {
185  collections.Append(new H1_FECollection((int)std::pow(2, level+1), dim));
186  fespaces.AddOrderRefinedLevel(collections.Last());
187  }
188 
189  cout << "Number of finite element unknowns: "
190  << fespaces.GetFinestFESpace().GetTrueVSize() << endl;
191 
192  // 6. Set up the linear form b(.) which corresponds to the right-hand side of
193  // the FEM linear system, which in this case is (1,phi_i) where phi_i are
194  // the basis functions in the finite element fespace.
195  LinearForm *b = new LinearForm(&fespaces.GetFinestFESpace());
196  ConstantCoefficient one(1.0);
197  b->AddDomainIntegrator(new DomainLFIntegrator(one));
198  b->Assemble();
199 
200  // 7. Define the solution vector x as a finite element grid function
201  // corresponding to fespace. Initialize x with initial guess of zero,
202  // which satisfies the boundary conditions.
203  GridFunction x(&fespaces.GetFinestFESpace());
204  x = 0.0;
205 
206  // 8. Create the multigrid operator using the previously created
207  // FiniteElementSpaceHierarchy and additional boundary information. This
208  // operator is then used to create the MultigridSolver as a preconditioner
209  // in the iterative solver.
210  Array<int> ess_bdr(mesh->bdr_attributes.Max());
211  ess_bdr = 1;
212 
213  DiffusionMultigrid M(fespaces, ess_bdr);
214  M.SetCycleType(Multigrid::CycleType::VCYCLE, 1, 1);
215 
216  OperatorPtr A;
217  Vector B, X;
218  M.FormFineLinearSystem(x, *b, A, X, B);
219  cout << "Size of linear system: " << A->Height() << endl;
220 
221  // 9. Solve the linear system A X = B.
222  PCG(*A, M, B, X, 1, 2000, 1e-12, 0.0);
223 
224  // 10. Recover the solution as a finite element grid function.
225  M.RecoverFineFEMSolution(X, *b, x);
226 
227  // 11. Save the refined mesh and the solution. This output can be viewed
228  // later using GLVis: "glvis -m refined.mesh -g sol.gf".
229  ofstream mesh_ofs("refined.mesh");
230  mesh_ofs.precision(8);
231  fespaces.GetFinestFESpace().GetMesh()->Print(mesh_ofs);
232  ofstream sol_ofs("sol.gf");
233  sol_ofs.precision(8);
234  x.Save(sol_ofs);
235 
236  // 12. Send the solution by socket to a GLVis server.
237  if (visualization)
238  {
239  char vishost[] = "localhost";
240  int visport = 19916;
241  socketstream sol_sock(vishost, visport);
242  sol_sock.precision(8);
243  sol_sock << "solution\n" << *fespaces.GetFinestFESpace().GetMesh() << x <<
244  flush;
245  }
246 
247  // 13. Free the used memory.
248  delete b;
249  for (int level = 0; level < collections.Size(); ++level)
250  {
251  delete collections[level];
252  }
253 
254  return 0;
255 }
virtual void AddOrderRefinedLevel(FiniteElementCollection *fec, int dim=1, int ordering=Ordering::byVDIM)
Adds one level to the hierarchy by using a different finite element order defined through FiniteEleme...
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:108
int visport
Conjugate gradient method.
Definition: solvers.hpp:493
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
Chebyshev accelerated smoothing with given vector, no matrix necessary.
Definition: solvers.hpp:379
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:84
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:331
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1020
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:462
Geometric multigrid associated with a hierarchy of finite element spaces.
Definition: multigrid.hpp:165
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:68
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Get a list of essential true dofs, ess_tdof_list, corresponding to the boundary attributes marked in ...
Definition: fespace.cpp:587
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Definition: device.cpp:279
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level.
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:159
STL namespace.
virtual const FiniteElementSpace & GetFinestFESpace() const
Returns the finite element space at the finest level.
int GetNumLevels() const
Returns the number of levels in the hierarchy.
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition: solvers.cpp:71
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
char vishost[]
int Append(const T &el)
Append element &#39;el&#39; to array, resize if necessary.
Definition: array.hpp:759
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:10232
void SetMaxIter(int max_it)
Definition: solvers.hpp:201
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:913
int main(int argc, char *argv[])
Definition: ex26.cpp:115
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:712
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:555
void SetAbsTol(double atol)
Definition: solvers.hpp:200
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:275
void SetRelTol(double rtol)
Definition: solvers.hpp:199
virtual const FiniteElementSpace & GetFESpaceAtLevel(int level) const
Returns the finite element space at the given level.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:219
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
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:1086
A "square matrix" operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
virtual void AddUniformlyRefinedLevel(int dim=1, int ordering=Ordering::byVDIM)
Adds one level to the hierarchy by uniformly refining the mesh on the previous level.
void SetOperatorOwner(bool own=true)
Set the ownership flag for the held Operator.
Definition: handle.hpp:120
int dim
Definition: ex24.cpp:53
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
Operator * Ptr() const
Access the underlying Operator pointer.
Definition: handle.hpp:87
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
T & Last()
Return the last element in the array.
Definition: array.hpp:792
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:507
Vector data type.
Definition: vector.hpp:58
virtual void Print(std::ostream &os=mfem::out) const
Definition: mesh.hpp:2011
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:259
Vector with associated FE space and LinearFormIntegrators.
Definition: linearform.hpp:24
Base class for solvers.
Definition: operator.hpp:682
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Definition: device.hpp:121
void SetType(Operator::Type tid)
Invoke Clear() and set a new type id.
Definition: handle.hpp:132