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