MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
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
38using namespace std;
39using namespace mfem;
40
41int 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 real_t 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
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 real_t 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}
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
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 if the AssemblyLevel is AssemblyLevel::LEGACY....
Adds new Domain Integrator. Assumes ownership of bfi.
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
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.
const SparseMatrix & SpMat() const
Returns a const reference to the sparse matrix: .
A class to handle Block diagonal preconditioners in a matrix-free implementation.
void SetDiagonalBlock(int iblock, Operator *op)
Add a square block op in the block-entry (iblock, iblock).
A class to handle Block systems in a matrix-free implementation.
void SetBlock(int iRow, int iCol, Operator *op, real_t c=1.0)
Add a block op in the block-entry (iblock, jblock).
virtual void Mult(const Vector &x, Vector &y) const
Operator application.
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator.
A class to handle Vectors in a block fashion.
Definition solvers.hpp:513
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition solvers.hpp:526
A coefficient that is constant across space and time.
Class for domain integration .
Definition lininteg.hpp:109
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition fe_coll.hpp:27
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:220
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition fespace.hpp:713
Class for grid function - Vector with associated FE space.
Definition gridfunc.hpp:31
virtual void MakeRef(FiniteElementSpace *f, real_t *v)
Make the GridFunction reference external data on a new FiniteElementSpace.
Definition gridfunc.cpp:203
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
Integrator that inverts the matrix assembled by another integrator.
void SetRelTol(real_t rtol)
Definition solvers.hpp:209
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition solvers.cpp:71
void SetMaxIter(int max_it)
Definition solvers.hpp:211
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition fe_coll.hpp:330
Vector with associated FE space and LinearFormIntegrators.
Adds new Domain Integrator. Assumes ownership of lfi.
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Mesh data type.
Definition mesh.hpp:56
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition mesh.hpp:282
virtual void Print(std::ostream &os=mfem::out, const std::string &comments="") const
Definition mesh.hpp:2288
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1226
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
int MeshGenerator() const
Get the mesh generator/type.
Definition mesh.hpp:1187
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization if the AssemblyLevel is AssemblyLevel::LEGACY.
void Assemble(int skip_zeros=1)
void EliminateTrialDofs(const Array< int > &bdr_attr_is_ess, const Vector &sol, Vector &rhs)
Eliminate essential boundary DOFs from the columns of the system.
const SparseMatrix & SpMat() const
Returns a const reference to the sparse matrix: .
Add a trace face integrator. Assumes ownership of bfi.
Adds a domain integrator. Assumes ownership of bfi.
Abstract operator.
Definition operator.hpp:25
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void PrintUsage(std::ostream &out) const
Print the usage message.
void PrintOptions(std::ostream &out) const
Print the options.
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 'var' to receive the value. Enable/disable tags are used to set the bool...
Definition optparser.hpp:82
bool Good() const
Return true if the command line options were parsed successfully.
The operator x -> R*A*P*x constructed through the actions of R^T, A and P.
Definition operator.hpp:817
Arbitrary order "H^{-1/2}-conforming" face finite elements defined on the interface between mesh elem...
Definition fe_coll.hpp:445
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition operator.hpp:686
Data type sparse matrix.
Definition sparsemat.hpp:51
real_t InnerProduct(const Vector &x, const Vector &y) const
Compute y^t A x.
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Integrator defining a sum of multiple Integrators.
Direct sparse solver using UMFPACK.
Definition solvers.hpp:1096
Vector data type.
Definition vector.hpp:80
int dim
Definition ex24.cpp:53
int main()
real_t b
Definition lissajous.cpp:42
const int visport
void RAP(const DenseMatrix &A, const DenseMatrix &P, DenseMatrix &RAP)
void PCG(const Operator &A, Solver &B, const Vector &b, Vector &x, int print_iter, int max_num_iter, real_t RTOLERANCE, real_t ATOLERANCE)
Preconditioned conjugate gradient method. (tolerances are squared)
Definition solvers.cpp:913
float real_t
Definition config.hpp:43
const char vishost[]