MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
ex5p.cpp
Go to the documentation of this file.
1// MFEM Example 5 - Parallel Version
2//
3// Compile with: make ex5p
4//
5// Sample runs: mpirun -np 4 ex5p -m ../data/square-disc.mesh
6// mpirun -np 4 ex5p -m ../data/star.mesh
7// mpirun -np 4 ex5p -m ../data/star.mesh -r 2 -pa
8// mpirun -np 4 ex5p -m ../data/beam-tet.mesh
9// mpirun -np 4 ex5p -m ../data/beam-hex.mesh
10// mpirun -np 4 ex5p -m ../data/beam-hex.mesh -pa
11// mpirun -np 4 ex5p -m ../data/escher.mesh
12// mpirun -np 4 ex5p -m ../data/fichera.mesh
13//
14// Device sample runs:
15// mpirun -np 4 ex5p -m ../data/star.mesh -r 2 -pa -d cuda
16// mpirun -np 4 ex5p -m ../data/star.mesh -r 2 -pa -d raja-cuda
17// mpirun -np 4 ex5p -m ../data/star.mesh -r 2 -pa -d raja-omp
18// mpirun -np 4 ex5p -m ../data/beam-hex.mesh -pa -d cuda
19//
20// Description: This example code solves a simple 2D/3D mixed Darcy problem
21// corresponding to the saddle point system
22//
23// k*u + grad p = f
24// - div u = g
25//
26// with natural boundary condition -p = <given pressure>.
27// Here, we use a given exact solution (u,p) and compute the
28// corresponding r.h.s. (f,g). We discretize with Raviart-Thomas
29// finite elements (velocity u) and piecewise discontinuous
30// polynomials (pressure p).
31//
32// The example demonstrates the use of the BlockOperator class, as
33// well as the collective saving of several grid functions in
34// VisIt (visit.llnl.gov) and ParaView (paraview.org) formats.
36// also illustrated.
37//
38// We recommend viewing examples 1-4 before viewing this example.
39
40#include "mfem.hpp"
41#include <fstream>
42#include <iostream>
43
44using namespace std;
45using namespace mfem;
46
47// Define the analytical solution and forcing terms / boundary conditions
48void uFun_ex(const Vector & x, Vector & u);
49real_t pFun_ex(const Vector & x);
50void fFun(const Vector & x, Vector & f);
51real_t gFun(const Vector & x);
52real_t f_natural(const Vector & x);
53
54int main(int argc, char *argv[])
55{
56 StopWatch chrono;
57
58 // 1. Initialize MPI and HYPRE.
59 Mpi::Init(argc, argv);
60 int num_procs = Mpi::WorldSize();
61 int myid = Mpi::WorldRank();
63 bool verbose = (myid == 0);
64
65 // 2. Parse command-line options.
66 const char *mesh_file = "../data/star.mesh";
67 int ref_levels = -1;
68 int order = 1;
69 bool par_format = false;
70 bool pa = false;
71 const char *device_config = "cpu";
72 bool visualization = 1;
74
75 OptionsParser args(argc, argv);
77 "Mesh file to use.");
79 "Number of times to refine the mesh uniformly.");
81 "Finite element order (polynomial degree).");
83 "--serial-format",
84 "Format to use when saving the results for VisIt.");
86 "--no-partial-assembly", "Enable Partial Assembly.");
88 "Device configuration string, see Device::Configure().");
90 "--no-visualization",
91 "Enable or disable GLVis visualization.");
94 "Save data using adios2 streams.");
95 args.Parse();
96 if (!args.Good())
97 {
98 if (verbose)
99 {
100 args.PrintUsage(cout);
101 }
102 return 1;
103 }
104 if (verbose)
105 {
106 args.PrintOptions(cout);
107 }
108
109 // 3. Enable hardware devices such as GPUs, and programming models such as
110 // CUDA, OCCA, RAJA and OpenMP based on command line options.
111 Device device(device_config);
112 if (myid == 0) { device.Print(); }
113
114 // 4. Read the (serial) mesh from the given mesh file on all processors. We
115 // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
116 // and volume meshes with the same code.
117 Mesh *mesh = new Mesh(mesh_file, 1, 1);
118 int dim = mesh->Dimension();
119
120 // 5. Refine the serial mesh on all processors to increase the resolution. In
121 // this example we do 'ref_levels' of uniform refinement. We choose
122 // 'ref_levels' to be the largest number that gives a final mesh with no
123 // more than 10,000 elements, unless the user specifies it as input.
124 {
125 if (ref_levels == -1)
126 {
127 ref_levels = (int)floor(log(10000./mesh->GetNE())/log(2.)/dim);
128 }
129
130 for (int l = 0; l < ref_levels; l++)
131 {
132 mesh->UniformRefinement();
133 }
134 }
135
136 // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
137 // this mesh further in parallel to increase the resolution. Once the
138 // parallel mesh is defined, the serial mesh can be deleted.
139 ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
140 delete mesh;
141 {
142 int par_ref_levels = 2;
143 for (int l = 0; l < par_ref_levels; l++)
144 {
145 pmesh->UniformRefinement();
146 }
147 }
148
149 // 7. Define a parallel finite element space on the parallel mesh. Here we
150 // use the Raviart-Thomas finite elements of the specified order.
151 FiniteElementCollection *hdiv_coll(new RT_FECollection(order, dim));
152 FiniteElementCollection *l2_coll(new L2_FECollection(order, dim));
153
154 ParFiniteElementSpace *R_space = new ParFiniteElementSpace(pmesh, hdiv_coll);
155 ParFiniteElementSpace *W_space = new ParFiniteElementSpace(pmesh, l2_coll);
156
157 HYPRE_BigInt dimR = R_space->GlobalTrueVSize();
158 HYPRE_BigInt dimW = W_space->GlobalTrueVSize();
159
160 if (verbose)
161 {
162 std::cout << "***********************************************************\n";
163 std::cout << "dim(R) = " << dimR << "\n";
164 std::cout << "dim(W) = " << dimW << "\n";
165 std::cout << "dim(R+W) = " << dimR + dimW << "\n";
166 std::cout << "***********************************************************\n";
167 }
168
169 // 8. Define the two BlockStructure of the problem. block_offsets is used
170 // for Vector based on dof (like ParGridFunction or ParLinearForm),
171 // block_trueOffstes is used for Vector based on trueDof (HypreParVector
172 // for the rhs and solution of the linear system). The offsets computed
173 // here are local to the processor.
174 Array<int> block_offsets(3); // number of variables + 1
175 block_offsets[0] = 0;
176 block_offsets[1] = R_space->GetVSize();
177 block_offsets[2] = W_space->GetVSize();
178 block_offsets.PartialSum();
179
180 Array<int> block_trueOffsets(3); // number of variables + 1
181 block_trueOffsets[0] = 0;
182 block_trueOffsets[1] = R_space->TrueVSize();
183 block_trueOffsets[2] = W_space->TrueVSize();
184 block_trueOffsets.PartialSum();
185
186 // 9. Define the coefficients, analytical solution, and rhs of the PDE.
187 ConstantCoefficient k(1.0);
188
192
195
196 // 10. Define the parallel grid function and parallel linear forms, solution
197 // vector and rhs.
198 MemoryType mt = device.GetMemoryType();
199 BlockVector x(block_offsets, mt), rhs(block_offsets, mt);
200 BlockVector trueX(block_trueOffsets, mt), trueRhs(block_trueOffsets, mt);
201
202 ParLinearForm *fform(new ParLinearForm);
203 fform->Update(R_space, rhs.GetBlock(0), 0);
206 fform->Assemble();
207 fform->SyncAliasMemory(rhs);
208 fform->ParallelAssemble(trueRhs.GetBlock(0));
209 trueRhs.GetBlock(0).SyncAliasMemory(trueRhs);
210
211 ParLinearForm *gform(new ParLinearForm);
212 gform->Update(W_space, rhs.GetBlock(1), 0);
214 gform->Assemble();
215 gform->SyncAliasMemory(rhs);
216 gform->ParallelAssemble(trueRhs.GetBlock(1));
217 trueRhs.GetBlock(1).SyncAliasMemory(trueRhs);
218
219 // 11. Assemble the finite element matrices for the Darcy operator
220 //
221 // D = [ M B^T ]
222 // [ B 0 ]
223 // where:
224 //
225 // M = \int_\Omega k u_h \cdot v_h d\Omega u_h, v_h \in R_h
226 // B = -\int_\Omega \div u_h q_h d\Omega u_h \in R_h, q_h \in W_h
227 ParBilinearForm *mVarf(new ParBilinearForm(R_space));
228 ParMixedBilinearForm *bVarf(new ParMixedBilinearForm(R_space, W_space));
229
230 HypreParMatrix *M = NULL;
231 HypreParMatrix *B = NULL;
232
233 if (pa) { mVarf->SetAssemblyLevel(AssemblyLevel::PARTIAL); }
235 mVarf->Assemble();
236 if (!pa) { mVarf->Finalize(); }
237
238 if (pa) { bVarf->SetAssemblyLevel(AssemblyLevel::PARTIAL); }
240 bVarf->Assemble();
241 if (!pa) { bVarf->Finalize(); }
242
243 BlockOperator *darcyOp = new BlockOperator(block_trueOffsets);
244
245 Array<int> empty_tdof_list; // empty
246 OperatorPtr opM, opB;
247
248 TransposeOperator *Bt = NULL;
249
250 if (pa)
251 {
252 mVarf->FormSystemMatrix(empty_tdof_list, opM);
253 bVarf->FormRectangularSystemMatrix(empty_tdof_list, empty_tdof_list, opB);
254 Bt = new TransposeOperator(opB.Ptr());
255
256 darcyOp->SetBlock(0,0, opM.Ptr());
257 darcyOp->SetBlock(0,1, Bt, -1.0);
258 darcyOp->SetBlock(1,0, opB.Ptr(), -1.0);
259 }
260 else
261 {
262 M = mVarf->ParallelAssemble();
263 B = bVarf->ParallelAssemble();
264 (*B) *= -1;
265 Bt = new TransposeOperator(B);
266
267 darcyOp->SetBlock(0,0, M);
268 darcyOp->SetBlock(0,1, Bt);
269 darcyOp->SetBlock(1,0, B);
270 }
271
272 // 12. Construct the operators for preconditioner
273 //
274 // P = [ diag(M) 0 ]
275 // [ 0 B diag(M)^-1 B^T ]
276 //
277 // Here we use Symmetric Gauss-Seidel to approximate the inverse of the
278 // pressure Schur Complement.
279 HypreParMatrix *MinvBt = NULL;
280 HypreParVector *Md = NULL;
281 HypreParMatrix *S = NULL;
282 Vector Md_PA;
283 Solver *invM, *invS;
284
285 if (pa)
286 {
287 Md_PA.SetSize(R_space->GetTrueVSize());
288 mVarf->AssembleDiagonal(Md_PA);
290 Vector invMd(Md_PA.Size());
291 for (int i=0; i<Md_PA.Size(); ++i)
292 {
293 invMd(i) = 1.0 / Md_host[i];
294 }
295
296 Vector BMBt_diag(W_space->GetTrueVSize());
298
299 Array<int> ess_tdof_list; // empty
300
301 invM = new OperatorJacobiSmoother(Md_PA, ess_tdof_list);
302 invS = new OperatorJacobiSmoother(BMBt_diag, ess_tdof_list);
303 }
304 else
305 {
306 Md = new HypreParVector(MPI_COMM_WORLD, M->GetGlobalNumRows(),
307 M->GetRowStarts());
308 M->GetDiag(*Md);
309
310 MinvBt = B->Transpose();
311 MinvBt->InvScaleRows(*Md);
312 S = ParMult(B, MinvBt);
313
314 invM = new HypreDiagScale(*M);
315 invS = new HypreBoomerAMG(*S);
316 }
317
318 invM->iterative_mode = false;
319 invS->iterative_mode = false;
320
322 block_trueOffsets);
323 darcyPr->SetDiagonalBlock(0, invM);
324 darcyPr->SetDiagonalBlock(1, invS);
325
326 // 13. Solve the linear system with MINRES.
327 // Check the norm of the unpreconditioned residual.
328 int maxIter(pa ? 1000 : 500);
329 real_t rtol(1.e-6);
330 real_t atol(1.e-10);
331
332 chrono.Clear();
333 chrono.Start();
334 MINRESSolver solver(MPI_COMM_WORLD);
335 solver.SetAbsTol(atol);
336 solver.SetRelTol(rtol);
337 solver.SetMaxIter(maxIter);
338 solver.SetOperator(*darcyOp);
339 solver.SetPreconditioner(*darcyPr);
340 solver.SetPrintLevel(verbose);
341 trueX = 0.0;
342 solver.Mult(trueRhs, trueX);
343 if (device.IsEnabled()) { trueX.HostRead(); }
344 chrono.Stop();
345
346 if (verbose)
347 {
348 if (solver.GetConverged())
349 std::cout << "MINRES converged in " << solver.GetNumIterations()
350 << " iterations with a residual norm of " << solver.GetFinalNorm() << ".\n";
351 else
352 std::cout << "MINRES did not converge in " << solver.GetNumIterations()
353 << " iterations. Residual norm is " << solver.GetFinalNorm() << ".\n";
354 std::cout << "MINRES solver took " << chrono.RealTime() << "s. \n";
355 }
356
357 // 14. Extract the parallel grid function corresponding to the finite element
358 // approximation X. This is the local solution on each processor. Compute
359 // L2 error norms.
362 u->MakeRef(R_space, x.GetBlock(0), 0);
363 p->MakeRef(W_space, x.GetBlock(1), 0);
364 u->Distribute(&(trueX.GetBlock(0)));
365 p->Distribute(&(trueX.GetBlock(1)));
366
367 int order_quad = max(2, 2*order+1);
369 for (int i=0; i < Geometry::NumGeom; ++i)
370 {
372 }
373
374 real_t err_u = u->ComputeL2Error(ucoeff, irs);
375 real_t norm_u = ComputeGlobalLpNorm(2, ucoeff, *pmesh, irs);
376 real_t err_p = p->ComputeL2Error(pcoeff, irs);
377 real_t norm_p = ComputeGlobalLpNorm(2, pcoeff, *pmesh, irs);
378
379 if (verbose)
380 {
381 std::cout << "|| u_h - u_ex || / || u_ex || = " << err_u / norm_u << "\n";
382 std::cout << "|| p_h - p_ex || / || p_ex || = " << err_p / norm_p << "\n";
383 }
384
385 // 15. Save the refined mesh and the solution in parallel. This output can be
386 // viewed later using GLVis: "glvis -np <np> -m mesh -g sol_*".
387 {
388 ostringstream mesh_name, u_name, p_name;
389 mesh_name << "mesh." << setfill('0') << setw(6) << myid;
390 u_name << "sol_u." << setfill('0') << setw(6) << myid;
391 p_name << "sol_p." << setfill('0') << setw(6) << myid;
392
393 ofstream mesh_ofs(mesh_name.str().c_str());
394 mesh_ofs.precision(8);
395 pmesh->Print(mesh_ofs);
396
397 ofstream u_ofs(u_name.str().c_str());
398 u_ofs.precision(8);
399 u->Save(u_ofs);
400
401 ofstream p_ofs(p_name.str().c_str());
402 p_ofs.precision(8);
403 p->Save(p_ofs);
404 }
405
406 // 16. Save data in the VisIt format
407 VisItDataCollection visit_dc("Example5-Parallel", pmesh);
408 visit_dc.RegisterField("velocity", u);
409 visit_dc.RegisterField("pressure", p);
410 visit_dc.SetFormat(!par_format ?
413 visit_dc.Save();
414
415 // 17. Save data in the ParaView format
416 ParaViewDataCollection paraview_dc("Example5P", pmesh);
417 paraview_dc.SetPrefixPath("ParaView");
418 paraview_dc.SetLevelsOfDetail(order);
419 paraview_dc.SetDataFormat(VTKFormat::BINARY);
420 paraview_dc.SetHighOrderOutput(true);
421 paraview_dc.SetCycle(0);
422 paraview_dc.SetTime(0.0);
423 paraview_dc.RegisterField("velocity",u);
424 paraview_dc.RegisterField("pressure",p);
425 paraview_dc.Save();
426
427 // 18. Optionally output a BP (binary pack) file using ADIOS2. This can be
428 // visualized with the ParaView VTX reader.
431 {
432 std::string postfix(mesh_file);
433 postfix.erase(0, std::string("../data/").size() );
434 postfix += "_o" + std::to_string(order);
435 const std::string collection_name = "ex5-p_" + postfix + ".bp";
436
444 }
445#endif
446
447 // 19. Send the solution by socket to a GLVis server.
448 if (visualization)
449 {
450 char vishost[] = "localhost";
451 int visport = 19916;
453 u_sock << "parallel " << num_procs << " " << myid << "\n";
454 u_sock.precision(8);
455 u_sock << "solution\n" << *pmesh << *u << "window_title 'Velocity'"
456 << endl;
457 // Make sure all ranks have sent their 'u' solution before initiating
458 // another set of GLVis connections (one from each rank):
459 MPI_Barrier(pmesh->GetComm());
461 p_sock << "parallel " << num_procs << " " << myid << "\n";
462 p_sock.precision(8);
463 p_sock << "solution\n" << *pmesh << *p << "window_title 'Pressure'"
464 << endl;
465 }
466
467 // 20. Free the used memory.
468 delete fform;
469 delete gform;
470 delete u;
471 delete p;
472 delete darcyOp;
473 delete darcyPr;
474 delete invM;
475 delete invS;
476 delete S;
477 delete Md;
478 delete MinvBt;
479 delete Bt;
480 delete B;
481 delete M;
482 delete mVarf;
483 delete bVarf;
484 delete W_space;
485 delete R_space;
486 delete l2_coll;
487 delete hdiv_coll;
488 delete pmesh;
489
490 return 0;
491}
492
493
494void uFun_ex(const Vector & x, Vector & u)
495{
496 real_t xi(x(0));
497 real_t yi(x(1));
498 real_t zi(0.0);
499 if (x.Size() == 3)
500 {
501 zi = x(2);
502 }
503
504 u(0) = - exp(xi)*sin(yi)*cos(zi);
505 u(1) = - exp(xi)*cos(yi)*cos(zi);
506
507 if (x.Size() == 3)
508 {
509 u(2) = exp(xi)*sin(yi)*sin(zi);
510 }
511}
512
513// Change if needed
515{
516 real_t xi(x(0));
517 real_t yi(x(1));
518 real_t zi(0.0);
519
520 if (x.Size() == 3)
521 {
522 zi = x(2);
523 }
524
525 return exp(xi)*sin(yi)*cos(zi);
526}
527
528void fFun(const Vector & x, Vector & f)
529{
530 f = 0.0;
531}
532
534{
535 if (x.Size() == 3)
536 {
537 return -pFun_ex(x);
538 }
539 else
540 {
541 return 0;
542 }
543}
544
546{
547 return (-pFun_ex(x));
548}
void SetLevelsOfDetail(const int levels_of_detail) noexcept
void PartialSum()
Fill the entries of the array with the cumulative sum of the entries.
Definition array.cpp:103
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level.
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.
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).
A class to handle Vectors in a block fashion.
Vector & GetBlock(int i)
Get the i-th vector in the block.
A coefficient that is constant across space and time.
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
void SetTime(real_t t)
Set physical time (for time-dependent simulations)
virtual void SetFormat(int fmt)
Set the desired output mesh and data format.
void SetPrefixPath(const std::string &prefix)
Set the path where the DataCollection will be saved.
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Definition device.hpp:123
static MemoryType GetMemoryType()
(DEPRECATED) Equivalent to GetDeviceMemoryType().
Definition device.hpp:278
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Definition device.cpp:286
static bool IsEnabled()
Return true if any backend other than Backend::CPU is enabled.
Definition device.hpp:247
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
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition fespace.hpp:713
A general function coefficient.
static const int NumGeom
Definition geom.hpp:42
The BoomerAMG solver in hypre.
Definition hypre.hpp:1691
Jacobi preconditioner in hypre.
Definition hypre.hpp:1481
Wrapper for hypre's ParCSR matrix class.
Definition hypre.hpp:388
void GetDiag(Vector &diag) const
Get the local diagonal of the matrix.
Definition hypre.cpp:1557
void InvScaleRows(const Vector &s)
Scale the local row i by 1./s(i)
Definition hypre.cpp:2135
HYPRE_BigInt * GetRowStarts() const
Return the parallel row partitioning array.
Definition hypre.hpp:689
HYPRE_BigInt GetGlobalNumRows() const
Return the global number of rows.
Definition hypre.hpp:679
HypreParMatrix * Transpose() const
Returns the transpose of *this.
Definition hypre.cpp:1684
Wrapper for hypre's parallel vector class.
Definition hypre.hpp:206
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
Definition hypre.hpp:74
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:100
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
real_t GetFinalNorm() const
Returns the final residual norm after termination of the solver during the last call to Mult().
Definition solvers.hpp:275
void SetRelTol(real_t rtol)
Definition solvers.hpp:209
int GetNumIterations() const
Returns the number of iterations taken during the last call to Mult()
Definition solvers.hpp:260
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
bool GetConverged() const
Returns true if the last call to Mult() converged successfully.
Definition solvers.hpp:262
void SetAbsTol(real_t atol)
Definition solvers.hpp:210
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition fe_coll.hpp:330
Adds new Domain Integrator. Assumes ownership of lfi.
Adds new Boundary Integrator. Assumes ownership of lfi.
MINRES method.
Definition solvers.hpp:628
virtual void Mult(const Vector &b, Vector &x) const
Iterative solution of the linear system using the MINRES method.
Definition solvers.cpp:1616
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition solvers.cpp:1595
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition solvers.hpp:640
Mesh data type.
Definition mesh.hpp:56
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
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 SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level. The default is AssemblyLevel::LEGACY.
void AssembleDiagonal_ADAt(const Vector &D, Vector &diag) const
Assemble the diagonal of ADA^T into diag, where A is this mixed bilinear form and D is a diagonal.
Adds a domain integrator. Assumes ownership of bfi.
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
static int WorldSize()
Return the size of MPI_COMM_WORLD.
static void Init(int &argc, char **&argv, int required=default_thread_required, int *provided=nullptr)
Singleton creation with Mpi::Init(argc, argv).
Pointer to an Operator of a specified type.
Definition handle.hpp:34
Operator * Ptr() const
Access the underlying Operator pointer.
Definition handle.hpp:87
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition solvers.hpp:313
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.
Class for parallel bilinear form.
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
void Assemble(int skip_zeros=1)
Assemble the local matrix.
virtual void AssembleDiagonal(Vector &diag) const
Assemble the diagonal of the bilinear form into diag. Note that diag is a true-dof Vector.
virtual void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)
Form the linear system matrix A, see FormLinearSystem() for details.
Abstract parallel finite element space.
Definition pfespace.hpp:29
HYPRE_BigInt GlobalTrueVSize() const
Definition pfespace.hpp:285
int GetTrueVSize() const override
Return the number of local vector true dofs.
Definition pfespace.hpp:289
int TrueVSize() const
Obsolete, kept for backward compatibility.
Definition pfespace.hpp:436
Class for parallel grid function.
Definition pgridfunc.hpp:33
Class for parallel linear form.
void Update(ParFiniteElementSpace *pf=NULL)
Update the object according to the given new FE space *pf.
void ParallelAssemble(Vector &tv)
Assemble the vector on the true dofs, i.e. P^t v.
void Assemble()
Assembles the ParLinearForm i.e. sums over all domain/bdr integrators.
Class for parallel meshes.
Definition pmesh.hpp:34
MPI_Comm GetComm() const
Definition pmesh.hpp:402
void Print(std::ostream &out=mfem::out, const std::string &comments="") const override
Definition pmesh.cpp:4801
Class for parallel bilinear form using different test and trial FE spaces.
virtual void FormRectangularSystemMatrix(const Array< int > &trial_tdof_list, const Array< int > &test_tdof_list, OperatorHandle &A)
Return in A a parallel (on truedofs) version of this operator.
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P_test^t A P_trial.
Helper class for ParaView visualization data.
void SetHighOrderOutput(bool high_order_output_)
void SetLevelsOfDetail(int levels_of_detail_)
virtual void Save() override
void SetDataFormat(VTKFormat fmt)
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition fe_coll.hpp:386
Base class for solvers.
Definition operator.hpp:683
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition operator.hpp:686
Timing object.
Definition tic_toc.hpp:36
double RealTime()
Return the number of real seconds elapsed since the stopwatch was started.
Definition tic_toc.cpp:432
void Start()
Start the stopwatch. The elapsed time is not cleared.
Definition tic_toc.cpp:411
void Stop()
Stop the stopwatch.
Definition tic_toc.cpp:422
void Clear()
Clear the elapsed time on the stopwatch and restart it if it's running.
Definition tic_toc.cpp:406
The transpose of a given operator. Switches the roles of the methods Mult() and MultTranspose().
Definition operator.hpp:750
for VectorFiniteElements (Nedelec, Raviart-Thomas)
Definition lininteg.hpp:347
A general vector function coefficient.
Vector data type.
Definition vector.hpp:80
virtual const real_t * HostRead() const
Definition vector.hpp:478
void SyncAliasMemory(const Vector &v) const
Update the alias memory location of the vector to match v.
Definition vector.hpp:259
int Size() const
Returns the size of the vector.
Definition vector.hpp:218
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:538
Data collection with VisIt I/O routines.
virtual void Save() override
Save the collection and a VisIt root file.
virtual void RegisterField(const std::string &field_name, GridFunction *gf) override
Add a grid function to the collection and update the root file.
int dim
Definition ex24.cpp:53
real_t pFun_ex(const Vector &x)
Definition ex5p.cpp:514
void uFun_ex(const Vector &x, Vector &u)
Definition ex5p.cpp:494
real_t gFun(const Vector &x)
Definition ex5p.cpp:533
real_t f_natural(const Vector &x)
Definition ex5p.cpp:545
void fFun(const Vector &x, Vector &f)
Definition ex5p.cpp:528
int main()
HYPRE_Int HYPRE_BigInt
const int visport
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
real_t ComputeGlobalLpNorm(real_t p, Coefficient &coeff, ParMesh &pmesh, const IntegrationRule *irs[])
Compute the global Lp norm of a function f. .
HypreParMatrix * ParMult(const HypreParMatrix *A, const HypreParMatrix *B, bool own_matrix)
Definition hypre.cpp:2949
float real_t
Definition config.hpp:43
MemoryType
Memory types supported by MFEM.
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30