MFEM  v4.3.0 Finite element discretization library
ex8p.cpp
Go to the documentation of this file.
1 // MFEM Example 8 - Parallel Version
2 //
3 // Compile with: make ex8p
4 //
5 // Sample runs: mpirun -np 4 ex8p -m ../data/square-disc.mesh
6 // mpirun -np 4 ex8p -m ../data/star.mesh
7 // mpirun -np 4 ex8p -m ../data/star-mixed.mesh
8 // mpirun -np 4 ex8p -m ../data/escher.mesh
9 // mpirun -np 4 ex8p -m ../data/fichera.mesh
10 // mpirun -np 4 ex8p -m ../data/fichera-mixed.mesh
11 // mpirun -np 4 ex8p -m ../data/square-disc-p2.vtk
12 // mpirun -np 4 ex8p -m ../data/square-disc-p3.mesh
13 // mpirun -np 4 ex8p -m ../data/star-surf.mesh -o 2
14 //
15 // Description: This example code demonstrates the use of the Discontinuous
16 // Petrov-Galerkin (DPG) method in its primal 2x2 block form as a
17 // simple finite element discretization of the Laplace problem
18 // -Delta u = f with homogeneous Dirichlet boundary conditions. We
19 // use high-order continuous trial space, a high-order interfacial
20 // (trace) space, and a high-order discontinuous test space
21 // defining a local dual (H^{-1}) norm.
22 //
23 // We use the primal form of DPG, see "A primal DPG method without
24 // a first-order reformulation", Demkowicz and Gopalakrishnan, CAM
25 // 2013, DOI:10.1016/j.camwa.2013.06.029.
26 //
27 // The example highlights the use of interfacial (trace) finite
28 // elements and spaces, trace face integrators and the definition
29 // of block operators and preconditioners. The use of the ADS
30 // preconditioner from hypre for interfacially-reduced H(div)
31 // problems is also illustrated.
32 //
33 // We recommend viewing examples 1-5 before viewing this example.
34
35 #include "mfem.hpp"
36 #include <fstream>
37 #include <iostream>
38
39 using namespace std;
40 using namespace mfem;
41
42 int main(int argc, char *argv[])
43 {
44  // 1. Initialize MPI.
45  int num_procs, myid;
46  MPI_Init(&argc, &argv);
47  MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
48  MPI_Comm_rank(MPI_COMM_WORLD, &myid);
49
50  // 2. Parse command-line options.
51  const char *mesh_file = "../data/star.mesh";
52  int order = 1;
53  bool visualization = 1;
54
55  OptionsParser args(argc, argv);
57  "Mesh file to use.");
59  "Finite element order (polynomial degree).");
61  "--no-visualization",
62  "Enable or disable GLVis visualization.");
63  args.Parse();
64  if (!args.Good())
65  {
66  if (myid == 0)
67  {
68  args.PrintUsage(cout);
69  }
70  MPI_Finalize();
71  return 1;
72  }
73  if (myid == 0)
74  {
75  args.PrintOptions(cout);
76  }
77
78  // 3. Read the (serial) mesh from the given mesh file on all processors. We
79  // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
80  // and volume meshes with the same code.
81  Mesh *mesh = new Mesh(mesh_file, 1, 1);
82  int dim = mesh->Dimension();
83
84  // 4. Refine the serial mesh on all processors to increase the resolution. In
85  // this example we do 'ref_levels' of uniform refinement. We choose
86  // 'ref_levels' to be the largest number that gives a final mesh with no
87  // more than 10,000 elements.
88  {
89  int ref_levels =
90  (int)floor(log(10000./mesh->GetNE())/log(2.)/dim);
91  for (int l = 0; l < ref_levels; l++)
92  {
93  mesh->UniformRefinement();
94  }
95  }
96
97  // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
98  // this mesh further in parallel to increase the resolution. Once the
99  // parallel mesh is defined, the serial mesh can be deleted.
100  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
101  delete mesh;
102  {
103  int par_ref_levels = 1;
104  for (int l = 0; l < par_ref_levels; l++)
105  {
106  pmesh->UniformRefinement();
107  }
108  }
109  pmesh->ReorientTetMesh();
110
111  // 6. Define the trial, interfacial (trace) and test DPG spaces:
112  // - The trial space, x0_space, contains the non-interfacial unknowns and
113  // has the essential BC.
114  // - The interfacial space, xhat_space, contains the interfacial unknowns
115  // and does not have essential BC.
116  // - The test space, test_space, is an enriched space where the enrichment
117  // degree may depend on the spatial dimension of the domain, the type of
118  // the mesh and the trial space order.
119  unsigned int trial_order = order;
120  unsigned int trace_order = order - 1;
121  unsigned int test_order = order; /* reduced order, full order is
122  (order + dim - 1) */
123  if (dim == 2 && (order%2 == 0 || (pmesh->MeshGenerator() & 2 && order > 1)))
124  {
125  test_order++;
126  }
127  if (test_order < trial_order)
128  {
129  if (myid == 0)
130  {
131  cerr << "Warning, test space not enriched enough to handle primal"
132  << " trial space\n";
133  }
134  }
135
136  FiniteElementCollection *x0_fec, *xhat_fec, *test_fec;
137
138  x0_fec = new H1_FECollection(trial_order, dim);
139  xhat_fec = new RT_Trace_FECollection(trace_order, dim);
140  test_fec = new L2_FECollection(test_order, dim);
141
142  ParFiniteElementSpace *x0_space, *xhat_space, *test_space;
143
144  x0_space = new ParFiniteElementSpace(pmesh, x0_fec);
145  xhat_space = new ParFiniteElementSpace(pmesh, xhat_fec);
146  test_space = new ParFiniteElementSpace(pmesh, test_fec);
147
148  HYPRE_BigInt glob_true_s0 = x0_space->GlobalTrueVSize();
149  HYPRE_BigInt glob_true_s1 = xhat_space->GlobalTrueVSize();
150  HYPRE_BigInt glob_true_s_test = test_space->GlobalTrueVSize();
151  if (myid == 0)
152  {
153  cout << "\nNumber of Unknowns:\n"
154  << " Trial space, X0 : " << glob_true_s0
155  << " (order " << trial_order << ")\n"
156  << " Interface space, Xhat : " << glob_true_s1
157  << " (order " << trace_order << ")\n"
158  << " Test space, Y : " << glob_true_s_test
159  << " (order " << test_order << ")\n\n";
160  }
161
162  // 7. Set up the linear form F(.) which corresponds to the right-hand side of
163  // the FEM linear system, which in this case is (f,phi_i) where f=1.0 and
164  // phi_i are the basis functions in the test finite element fespace.
165  ConstantCoefficient one(1.0);
166  ParLinearForm * F = new ParLinearForm(test_space);
168  F->Assemble();
169
170  ParGridFunction * x0 = new ParGridFunction(x0_space);
171  *x0 = 0.;
172
173  // 8. Set up the mixed bilinear form for the primal trial unknowns, B0,
174  // the mixed bilinear form for the interfacial unknowns, Bhat,
175  // the inverse stiffness matrix on the discontinuous test space, Sinv,
176  // and the stiffness matrix on the continuous trial space, S0.
177  Array<int> ess_bdr(pmesh->bdr_attributes.Max());
178  ess_bdr = 1;
179  Array<int> ess_dof;
180  x0_space->GetEssentialVDofs(ess_bdr, ess_dof);
181
182  ParMixedBilinearForm *B0 = new ParMixedBilinearForm(x0_space,test_space);
184  B0->Assemble();
185  B0->EliminateEssentialBCFromTrialDofs(ess_dof, *x0, *F);
186  B0->Finalize();
187
188  ParMixedBilinearForm *Bhat = new ParMixedBilinearForm(xhat_space,test_space);
190  Bhat->Assemble();
191  Bhat->Finalize();
192
193  ParBilinearForm *Sinv = new ParBilinearForm(test_space);
194  SumIntegrator *Sum = new SumIntegrator;
198  Sinv->Assemble();
199  Sinv->Finalize();
200
201  ParBilinearForm *S0 = new ParBilinearForm(x0_space);
203  S0->Assemble();
204  S0->EliminateEssentialBC(ess_bdr);
205  S0->Finalize();
206
207  HypreParMatrix * matB0 = B0->ParallelAssemble(); delete B0;
208  HypreParMatrix * matBhat = Bhat->ParallelAssemble(); delete Bhat;
209  HypreParMatrix * matSinv = Sinv->ParallelAssemble(); delete Sinv;
210  HypreParMatrix * matS0 = S0->ParallelAssemble(); delete S0;
211
212  // 9. Define the block structure of the problem, by creating the offset
213  // variables. Also allocate two BlockVector objects to store the solution
214  // and rhs.
215  enum {x0_var, xhat_var, NVAR};
216
217  int true_s0 = x0_space->TrueVSize();
218  int true_s1 = xhat_space->TrueVSize();
219  int true_s_test = test_space->TrueVSize();
220
221  Array<int> true_offsets(NVAR+1);
222  true_offsets[0] = 0;
223  true_offsets[1] = true_s0;
224  true_offsets[2] = true_s0+true_s1;
225
226  Array<int> true_offsets_test(2);
227  true_offsets_test[0] = 0;
228  true_offsets_test[1] = true_s_test;
229
230  BlockVector x(true_offsets), b(true_offsets);
231  x = 0.0;
232  b = 0.0;
233
234  // 10. Set up the 1x2 block Least Squares DPG operator, B = [B0 Bhat],
235  // the normal equation operator, A = B^t Sinv B, and
236  // the normal equation right-hand-size, b = B^t Sinv F.
237  BlockOperator B(true_offsets_test, true_offsets);
238  B.SetBlock(0, 0, matB0);
239  B.SetBlock(0, 1, matBhat);
240
241  RAPOperator A(B, *matSinv, B);
242
243  HypreParVector *trueF = F->ParallelAssemble();
244  {
245  HypreParVector SinvF(test_space);
246  matSinv->Mult(*trueF, SinvF);
247  B.MultTranspose(SinvF, b);
248  }
249
250  // 11. Set up a block-diagonal preconditioner for the 2x2 normal equation
251  //
252  // [ S0^{-1} 0 ]
253  // [ 0 Shat^{-1} ] Shat = (Bhat^T Sinv Bhat)
254  //
255  // corresponding to the primal (x0) and interfacial (xhat) unknowns.
256  // Since the Shat operator is equivalent to an H(div) matrix reduced to
257  // the interfacial skeleton, we approximate its inverse with one V-cycle
258  // of the ADS preconditioner from the hypre library (in 2D we use AMS for
259  // the rotated H(curl) problem).
260  HypreBoomerAMG *S0inv = new HypreBoomerAMG(*matS0);
261  S0inv->SetPrintLevel(0);
262
263  HypreParMatrix *Shat = RAP(matSinv, matBhat);
264  HypreSolver *Shatinv;
265  if (dim == 2) { Shatinv = new HypreAMS(*Shat, xhat_space); }
266  else { Shatinv = new HypreADS(*Shat, xhat_space); }
267
268  BlockDiagonalPreconditioner P(true_offsets);
269  P.SetDiagonalBlock(0, S0inv);
270  P.SetDiagonalBlock(1, Shatinv);
271
272  // 12. Solve the normal equation system using the PCG iterative solver.
273  // Check the weighted norm of residual for the DPG least square problem.
274  // Wrap the primal variable in a GridFunction for visualization purposes.
275  CGSolver pcg(MPI_COMM_WORLD);
276  pcg.SetOperator(A);
277  pcg.SetPreconditioner(P);
278  pcg.SetRelTol(1e-6);
279  pcg.SetMaxIter(200);
280  pcg.SetPrintLevel(1);
281  pcg.Mult(b, x);
282
283  {
284  HypreParVector LSres(test_space), tmp(test_space);
285  B.Mult(x, LSres);
286  LSres -= *trueF;
287  matSinv->Mult(LSres, tmp);
288  double res = sqrt(InnerProduct(LSres, tmp));
289  if (myid == 0)
290  {
291  cout << "\n|| B0*x0 + Bhat*xhat - F ||_{S^-1} = " << res << endl;
292  }
293  }
294
295  x0->Distribute(x.GetBlock(x0_var));
296
297  // 13. Save the refined mesh and the solution in parallel. This output can
298  // be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
299  {
300  ostringstream mesh_name, sol_name;
301  mesh_name << "mesh." << setfill('0') << setw(6) << myid;
302  sol_name << "sol." << setfill('0') << setw(6) << myid;
303
304  ofstream mesh_ofs(mesh_name.str().c_str());
305  mesh_ofs.precision(8);
306  pmesh->Print(mesh_ofs);
307
308  ofstream sol_ofs(sol_name.str().c_str());
309  sol_ofs.precision(8);
310  x0->Save(sol_ofs);
311  }
312
313  // 14. Send the solution by socket to a GLVis server.
314  if (visualization)
315  {
316  char vishost[] = "localhost";
317  int visport = 19916;
318  socketstream sol_sock(vishost, visport);
319  sol_sock << "parallel " << num_procs << " " << myid << "\n";
320  sol_sock.precision(8);
321  sol_sock << "solution\n" << *pmesh << *x0 << flush;
322  }
323
324  // 15. Free the used memory.
325  delete trueF;
326  delete Shatinv;
327  delete S0inv;
328  delete Shat;
329  delete matB0;
330  delete matBhat;
331  delete matSinv;
332  delete matS0;
333  delete x0;
334  delete F;
335  delete test_space;
336  delete xhat_space;
337  delete x0_space;
338  delete test_fec;
339  delete xhat_fec;
340  delete x0_fec;
341  delete pmesh;
342
343  MPI_Finalize();
344
345  return 0;
346 }
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:97
Definition: solvers.hpp:316
The Auxiliary-space Maxwell Solver in hypre.
Definition: hypre.hpp:1540
Integrator defining a sum of multiple Integrators.
Definition: bilininteg.hpp:363
The Auxiliary-space Divergence Solver in hypre.
Definition: hypre.hpp:1579
A class to handle Vectors in a block fashion.
Definition: blockvector.hpp:30
virtual void GetEssentialVDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_dofs, int component=-1) const
Determine the boundary degrees of freedom.
Definition: pfespace.cpp:776
A coefficient that is constant across space and time.
Definition: coefficient.hpp:78
HypreParMatrix * RAP(const HypreParMatrix *A, const HypreParMatrix *P)
Returns the matrix P^t * A * P.
Definition: hypre.cpp:2485
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:616
virtual void ReorientTetMesh()
See the remarks for the serial version in mesh.hpp.
Definition: pmesh.cpp:2798
int TrueVSize() const
Obsolete, kept for backward compatibility.
Definition: pfespace.hpp:422
HYPRE_BigInt GlobalTrueVSize() const
Definition: pfespace.hpp:275
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator.
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:846
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:841
Abstract parallel finite element space.
Definition: pfespace.hpp:28
Definition: bilininteg.hpp:375
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P_test^t A P_trial.
Adds a domain integrator. Assumes ownership of bfi.
Integrator that inverts the matrix assembled by another integrator.
Definition: bilininteg.hpp:343
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1387
Class for parallel linear form.
Definition: plinearform.hpp:26
void SetPrintLevel(int print_lvl)
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:150
A class to handle Block diagonal preconditioners in a matrix-free implementation. ...
constexpr char vishost[]
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9143
void SetPrintLevel(int print_level)
Definition: hypre.hpp:1467
constexpr int visport
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:100
int MeshGenerator()
Get the mesh generator/type.
Definition: mesh.hpp:839
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
void Assemble(int skip_zeros=1)
void Assemble(int skip_zeros=1)
Assemble the local matrix.
virtual void Print(std::ostream &out=mfem::out) const
Definition: pmesh.cpp:4382
The operator x -&gt; R*A*P*x constructed through the actions of R^T, A and P.
Definition: operator.hpp:756
int Dimension() const
Definition: mesh.hpp:911
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:457
Add a trace face integrator. Assumes ownership of bfi.
virtual void Mult(const Vector &x, Vector &y) const
Operator application.
Wrapper for hypre&#39;s parallel vector class.
Definition: hypre.hpp:99
Adds new Domain Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:39
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:204
void SetRelTol(double rtol)
Definition: solvers.hpp:98
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
HYPRE_Int HYPRE_BigInt
void Distribute(const Vector *tv)
Definition: pgridfunc.cpp:132
void EliminateEssentialBCFromTrialDofs(const Array< int > &marked_vdofs, const Vector &sol, Vector &rhs)
double InnerProduct(HypreParVector *x, HypreParVector *y)
Definition: hypre.cpp:316
Class for parallel bilinear form using different test and trial FE spaces.
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
int dim
Definition: ex24.cpp:53
Adds new Domain Integrator. Assumes ownership of bfi.
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:327
Class for parallel bilinear form.
Abstract class for hypre&#39;s solvers and preconditioners.
Definition: hypre.hpp:970
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:330
Arbitrary order &quot;H^{-1/2}-conforming&quot; face finite elements defined on the interface between mesh elem...
Definition: fe_coll.hpp:396
HYPRE_Int Mult(HypreParVector &x, HypreParVector &y, double alpha=1.0, double beta=0.0) const
Computes y = alpha * A * x + beta * y.
Definition: hypre.cpp:1529
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:92
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:216
Class for parallel grid function.
Definition: pgridfunc.hpp:32
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:277
A class to handle Block systems in a matrix-free implementation.
Class for parallel meshes.
Definition: pmesh.hpp:32
void ParallelAssemble(Vector &tv)
Assemble the vector on the true dofs, i.e. P^t v.
Definition: plinearform.cpp:87
int main()
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 &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:285
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150