MFEM  v3.0
ex8.cpp
Go to the documentation of this file.
1 // MFEM Example 8
2 //
3 // Compile with: make ex8
4 //
5 // Sample runs: ex8 -m ../data/square-disc.mesh
6 // ex8 -m ../data/star.mesh
7 // ex8 -m ../data/escher.mesh
8 // ex8 -m ../data/fichera.mesh
9 // ex8 -m ../data/square-disc-p2.vtk
10 // ex8 -m ../data/square-disc-p3.mesh
11 // ex8 -m ../data/star-surf.mesh -o 2
12 //
13 // Description: This example code demonstrates the use of the Discontinuous
14 // Petrov-Galerkin (DPG) method in its primal 2x2 block form as a
15 // simple finite element discretization of the Laplace problem
16 // -Delta u = f with homogeneous Dirichlet boundary conditions. We
17 // use high-order continuous trial space, a high-order interfacial
18 // (trace) space, and a high-order discontinuous test space
19 // defining a local dual (H^{-1}) norm.
20 //
21 // The example highlights the use of interfacial (trace) finite
22 // elements and spaces, trace face integrators and the definition
23 // of block operators and preconditioners.
24 //
25 // We recommend viewing examples 1-5 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 int main(int argc, char *argv[])
35 {
36  // 1. Parse command-line options.
37  const char *mesh_file = "../data/star.mesh";
38  int order = 1;
39  bool visualization = 1;
40
41  OptionsParser args(argc, argv);
43  "Mesh file to use.");
45  "Finite element order (polynomial degree).");
47  "--no-visualization",
48  "Enable or disable GLVis visualization.");
49  args.Parse();
50  if (!args.Good())
51  {
52  args.PrintUsage(cout);
53  return 1;
54  }
55  args.PrintOptions(cout);
56
57  // 2. Read the mesh from the given mesh file. We can handle triangular,
58  // quadrilateral, tetrahedral, hexahedral, surface and volume meshes with
59  // the same code.
60  Mesh *mesh;
61  ifstream imesh(mesh_file);
62  if (!imesh)
63  {
64  cerr << "\nCan not open mesh file: " << mesh_file << '\n' << endl;
65  return 2;
66  }
67  mesh = new Mesh(imesh, 1, 1);
68  imesh.close();
69  int dim = mesh->Dimension();
70
71  // 3. Refine the mesh to increase the resolution. In this example we do
72  // 'ref_levels' of uniform refinement. We choose 'ref_levels' to be the
73  // largest number that gives a final mesh with no more than 10,000
74  // elements.
75  {
76  int ref_levels =
77  (int)floor(log(10000./mesh->GetNE())/log(2.)/dim);
78  for (int l = 0; l < ref_levels; l++)
79  mesh->UniformRefinement();
80  }
81
82  // 4. Define the trial, interfacial (trace) and test DPG spaces:
83  // - The trial space, x0_space, contains the non-interfacial unknowns and
84  // has the essential BC.
85  // - The interfacial space, xhat_space, contains the interfacial unknowns
86  // and does not have essential BC.
87  // - The test space, test_space, is an enriched space where the enrichment
88  // degree may depend on the spatial dimension of the domain.
89  int trial_order = order;
90  int trace_order = order - 1;
91  int test_order = order; // reduced order, full order is (order + dim - 1)
92  if (test_order < trial_order)
93  cerr << "Warning, test space not enriched enough to handle primal"
94  << " trial space\n";
95
96  FiniteElementCollection *x0_fec, *xhat_fec, *test_fec;
97
98  x0_fec = new H1_FECollection(trial_order, dim);
99  xhat_fec = new RT_Trace_FECollection(trace_order, dim);
100  test_fec = new L2_FECollection(test_order, dim);
101
102  FiniteElementSpace *x0_space = new FiniteElementSpace(mesh, x0_fec);
103  FiniteElementSpace *xhat_space = new FiniteElementSpace(mesh, xhat_fec);
104  FiniteElementSpace *test_space = new FiniteElementSpace(mesh, test_fec);
105
106  // 5. Define the block structure of the problem, by creating the offset
107  // variables. Also allocate two BlockVector objects to store the solution
108  // and rhs.
109  enum {x0_var, xhat_var, NVAR};
110
111  int s0 = x0_space->GetVSize();
112  int s1 = xhat_space->GetVSize();
113  int s_test = test_space->GetVSize();
114
115  Array<int> offsets(NVAR+1);
116  offsets[0] = 0;
117  offsets[1] = s0;
118  offsets[2] = s0+s1;
119
120  Array<int> offsets_test(2);
121  offsets_test[0] = 0;
122  offsets_test[1] = s_test;
123
124  std::cout << "\nNumber of Unknowns:\n"
125  << " Trial space, X0 : " << s0 << '\n'
126  << " Interface space, Xhat : " << s1 << '\n'
127  << " Test space, Y : " << s_test << "\n\n";
128
129  BlockVector x(offsets), b(offsets);
130  x = 0.;
131
132  // 6. Set up the linear form F(.) which corresponds to the right-hand side of
133  // the FEM linear system, which in this case is (f,phi_i) where f=1.0 and
134  // phi_i are the basis functions in the test finite element fespace.
135  ConstantCoefficient one(1.0);
136  LinearForm F(test_space);
138  F.Assemble();
139
140  // 7. Set up the mixed bilinear form for the primal trial unknowns, B0,
141  // the mixed bilinear form for the interfacial unknowns, Bhat,
142  // the inverse stiffness matrix on the discontinuous test space, Sinv,
143  // and the stiffness matrix on the continuous trial space, S0.
144  Array<int> ess_bdr(mesh->bdr_attributes.Max());
145  ess_bdr = 1;
146
147  MixedBilinearForm *B0 = new MixedBilinearForm(x0_space,test_space);
149  B0->Assemble();
150  B0->EliminateTrialDofs(ess_bdr, x.GetBlock(x0_var), F);
151  B0->Finalize();
152
153  MixedBilinearForm *Bhat = new MixedBilinearForm(xhat_space,test_space);
155  Bhat->Assemble();
156  Bhat->Finalize();
157
158  BilinearForm *Sinv = new BilinearForm(test_space);
159  SumIntegrator *Sum = new SumIntegrator;
163  Sinv->Assemble();
164  Sinv->Finalize();
165
166  BilinearForm *S0 = new BilinearForm(x0_space);
169  S0->Assemble();
170  S0->EliminateEssentialBC(ess_bdr);
171  S0->Finalize();
172
173  SparseMatrix &matB0 = B0->SpMat();
174  SparseMatrix &matBhat = Bhat->SpMat();
175  SparseMatrix &matSinv = Sinv->SpMat();
176  SparseMatrix &matS0 = S0->SpMat();
177
178  // 8. Set up the 1x2 block Least Squares DPG operator, B = [B0 Bhat],
179  // the normal equation operator, A = B^t Sinv B, and
180  // the normal equation right-hand-size, b = B^t Sinv F.
181  BlockOperator B(offsets_test, offsets);
182  B.SetBlock(0,0,&matB0);
183  B.SetBlock(0,1,&matBhat);
184  RAPOperator A(B, matSinv, B);
185  {
186  Vector SinvF(s_test);
187  matSinv.Mult(F,SinvF);
188  B.MultTranspose(SinvF, b);
189  }
190
191  // 9. Set up a block-diagonal preconditioner for the 2x2 normal equation
192  //
193  // [ S0^{-1} 0 ]
194  // [ 0 Shat^{-1} ] Shat = (Bhat^T Sinv Bhat)
195  //
196  // corresponding to the primal (x0) and interfacial (xhat) unknowns.
197  SparseMatrix * Shat = RAP(matBhat, matSinv, matBhat);
198
199 #ifndef MFEM_USE_SUITESPARSE
200  const double prec_rtol = 1e-3;
201  const int prec_maxit = 200;
202  CGSolver *S0inv = new CGSolver;
203  S0inv->SetOperator(matS0);
204  S0inv->SetPrintLevel(-1);
205  S0inv->SetRelTol(prec_rtol);
206  S0inv->SetMaxIter(prec_maxit);
207  CGSolver *Shatinv = new CGSolver;
208  Shatinv->SetOperator(*Shat);
209  Shatinv->SetPrintLevel(-1);
210  Shatinv->SetRelTol(prec_rtol);
211  Shatinv->SetMaxIter(prec_maxit);
212  // Disable 'iterative_mode' when using CGSolver (or any IterativeSolver) as
213  // a preconditioner:
214  S0inv->iterative_mode = false;
215  Shatinv->iterative_mode = false;
216 #else
217  Operator *S0inv = new UMFPackSolver(matS0);
218  Operator *Shatinv = new UMFPackSolver(*Shat);
219 #endif
220
221  BlockDiagonalPreconditioner P(offsets);
222  P.SetDiagonalBlock(0, S0inv);
223  P.SetDiagonalBlock(1, Shatinv);
224
225  // 10. Solve the normal equation system using the PCG iterative solver.
226  // Check the weighted norm of residual for the DPG least square problem.
227  // Wrap the primal variable in a GridFunction for visualization purposes.
228  PCG(A, P, b, x, 1, 200, 1e-12, 0.0);
229
230  {
231  Vector LSres(s_test);
232  B.Mult(x, LSres);
233  LSres -= F;
234  double res = sqrt(matSinv.InnerProduct(LSres, LSres));
235  cout << "\n|| B0*x0 + Bhat*xhat - F ||_{S^-1} = " << res << endl;
236  }
237
238  GridFunction x0;
239  x0.Update(x0_space, x.GetBlock(x0_var), 0);
240
241  // 11. Save the refined mesh and the solution. This output can be viewed
242  // later using GLVis: "glvis -m refined.mesh -g sol.gf".
243  {
244  ofstream mesh_ofs("refined.mesh");
245  mesh_ofs.precision(8);
246  mesh->Print(mesh_ofs);
247  ofstream sol_ofs("sol.gf");
248  sol_ofs.precision(8);
249  x0.Save(sol_ofs);
250  }
251
252  // 12. Send the solution by socket to a GLVis server.
253  if (visualization)
254  {
255  char vishost[] = "localhost";
256  int visport = 19916;
257  socketstream sol_sock(vishost, visport);
258  sol_sock.precision(8);
259  sol_sock << "solution\n" << *mesh << x0 << flush;
260  }
261
262  // 13. Free the used memory.
263  delete S0inv;
264  delete Shatinv;
265  delete Shat;
266  delete Bhat;
267  delete B0;
268  delete S0;
269  delete Sinv;
270  delete test_space;
271  delete test_fec;
272  delete xhat_space;
273  delete xhat_fec;
274  delete x0_space;
275  delete x0_fec;
276  delete mesh;
277
278  return 0;
279 }
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:47
int GetVSize() const
Definition: fespace.hpp:151
Definition: solvers.hpp:109
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:26
Integrator defining a sum of multiple Integrators.
Definition: bilininteg.hpp:139
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
Subclass constant coefficient.
Definition: coefficient.hpp:57
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:34
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator.
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:396
Definition: bilininteg.hpp:149
void SetBlock(int iRow, int iCol, Operator *op)
Add a block op in the block-entry (iblock, jblock).
bool iterative_mode
If true, use the second argument of Mult as an initial guess.
Definition: operator.hpp:106
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
Direct sparse solver using UMFPACK.
Definition: solvers.hpp:329
Integrator that inverts the matrix assembled by another integrator.
Definition: bilininteg.hpp:121
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:1797
void SetPrintLevel(int print_lvl)
Definition: solvers.cpp:71
A class to handle Block diagonal preconditioners in a matrix-free implementation. ...
Data type sparse matrix.
Definition: sparsemat.hpp:38
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:6225
void SetMaxIter(int max_it)
Definition: solvers.hpp:61
HypreParMatrix * RAP(HypreParMatrix *A, HypreParMatrix *P)
Returns the matrix P^t * A * P.
Definition: hypre.cpp:781
T Max() const
Definition: array.cpp:78
void Assemble(int skip_zeros=1)
The operator x -&gt; R*A*P*x.
Definition: operator.hpp:164
void EliminateTrialDofs(Array< int > &bdr_attr_is_ess, Vector &sol, Vector &rhs)
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
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:385
virtual void Mult(const Vector &x, Vector &y) const
Operator application.
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
void SetRelTol(double rtol)
Definition: solvers.hpp:59
Abstract finite element space.
Definition: fespace.hpp:61
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Definition: sparsemat.cpp:312
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
const SparseMatrix & SpMat() const
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:266
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:123
Vector data type.
Definition: vector.hpp:29
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:52
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.
Abstract operator.
Definition: operator.hpp:21
A class to handle Block systems in a matrix-free implementation.
double InnerProduct(const Vector &x, const Vector &y) const
Compute y^t A x.
Definition: sparsemat.cpp:468
void EliminateEssentialBC(Array< int > &bdr_attr_is_ess, Vector &sol, Vector &rhs, int d=0)
void SetDiagonalBlock(int iblock, Operator *op)
Add a square block op in the block-entry (iblock, iblock).
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:83
bool Good() const
Definition: optparser.hpp:120