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// PETSc Modification
3//
4// Compile with: make ex5p
5//
6// Sample runs:
7// mpirun -np 4 ex5p -m ../../data/beam-tet.mesh --petscopts rc_ex5p_fieldsplit
8// mpirun -np 4 ex5p -m ../../data/star.mesh --petscopts rc_ex5p_bddc --nonoverlapping
9//
10// Description: This example code solves a simple 2D/3D mixed Darcy problem
11// corresponding to the saddle point system
12// k*u + grad p = f
13// - div u = g
14// with natural boundary condition -p = <given pressure>.
15// Here, we use a given exact solution (u,p) and compute the
16// corresponding r.h.s. (f,g). We discretize with Raviart-Thomas
17// finite elements (velocity u) and piecewise discontinuous
18// polynomials (pressure p).
19//
20// The example demonstrates the use of the BlockOperator class, as
21// well as the collective saving of several grid functions in a
22// VisIt (visit.llnl.gov) visualization format.
23//
24// Two types of PETSc solvers can be used: BDDC or fieldsplit.
25// When using BDDC, the nonoverlapping assembly feature should be
26// used. This specific example needs PETSc compiled with support
27// for SuiteSparse and/or MUMPS for using BDDC.
28//
29// We recommend viewing examples 1-4 before viewing this example.
30
31#include "mfem.hpp"
32#include <fstream>
33#include <iostream>
34
35#ifndef MFEM_USE_PETSC
36#error This example requires that MFEM is built with MFEM_USE_PETSC=YES
37#endif
38
39using namespace std;
40using namespace mfem;
41
42// Define the analytical solution and forcing terms / boundary conditions
43void uFun_ex(const Vector & x, Vector & u);
44real_t pFun_ex(const Vector & x);
45void fFun(const Vector & x, Vector & f);
46real_t gFun(const Vector & x);
47real_t f_natural(const Vector & x);
48
49int main(int argc, char *argv[])
50{
51 StopWatch chrono;
52
53 // 1. Initialize MPI and HYPRE.
54 Mpi::Init(argc, argv);
55 int num_procs = Mpi::WorldSize();
56 int myid = Mpi::WorldRank();
58 bool verbose = (myid == 0);
59
60 // 2. Parse command-line options.
61 const char *mesh_file = "../../data/star.mesh";
62 int order = 1;
63 bool par_format = false;
64 bool visualization = 1;
65 bool use_petsc = true;
66 bool use_nonoverlapping = false;
67 bool local_bdr_spec = false;
68 const char *petscrc_file = "";
69
70 OptionsParser args(argc, argv);
72 "Mesh file to use.");
74 "Finite element order (polynomial degree).");
76 "--serial-format",
77 "Format to use when saving the results for VisIt.");
79 "--no-visualization",
80 "Enable or disable GLVis visualization.");
82 "--no-petsc",
83 "Use or not PETSc to solve the linear system.");
85 "PetscOptions file to use.");
87 "-no-nonoverlapping", "--no-nonoverlapping",
88 "Use or not the block diagonal PETSc's matrix format "
89 "for non-overlapping domain decomposition.");
91 "--no-local-bdr",
92 "Specify boundary dofs in local (Vdofs) ordering.");
93 args.Parse();
94 if (!args.Good())
95 {
96 if (verbose)
97 {
98 args.PrintUsage(cout);
99 }
100 return 1;
101 }
102 if (verbose)
103 {
104 args.PrintOptions(cout);
105 }
106 // 2b. We initialize PETSc
107 if (use_petsc) { MFEMInitializePetsc(NULL,NULL,petscrc_file,NULL); }
108
109 // 3. Read the (serial) mesh from the given mesh file on all processors. We
110 // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
111 // and volume meshes with the same code.
112 Mesh *mesh = new Mesh(mesh_file, 1, 1);
113 int dim = mesh->Dimension();
114
115 // 4. Refine the serial mesh on all processors to increase the resolution. In
116 // this example we do 'ref_levels' of uniform refinement. We choose
117 // 'ref_levels' to be the largest number that gives a final mesh with no
118 // more than 10,000 elements.
119 {
120 int ref_levels =
121 (int)floor(log(10000./mesh->GetNE())/log(2.)/dim);
122 for (int l = 0; l < ref_levels; l++)
123 {
124 mesh->UniformRefinement();
125 }
126 }
127
128 // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
129 // this mesh further in parallel to increase the resolution. Once the
130 // parallel mesh is defined, the serial mesh can be deleted.
131 ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
132 delete mesh;
133 {
134 int par_ref_levels = 2;
135 for (int l = 0; l < par_ref_levels; l++)
136 {
137 pmesh->UniformRefinement();
138 }
139 }
140
141 // 6. Define a parallel finite element space on the parallel mesh. Here we
142 // use the Raviart-Thomas finite elements of the specified order.
143 FiniteElementCollection *hdiv_coll(new RT_FECollection(order, dim));
144 FiniteElementCollection *l2_coll(new L2_FECollection(order, dim));
145
146 ParFiniteElementSpace *R_space = new ParFiniteElementSpace(pmesh, hdiv_coll);
147 ParFiniteElementSpace *W_space = new ParFiniteElementSpace(pmesh, l2_coll);
148
149 HYPRE_BigInt dimR = R_space->GlobalTrueVSize();
150 HYPRE_BigInt dimW = W_space->GlobalTrueVSize();
151
152 if (verbose)
153 {
154 std::cout << "***********************************************************\n";
155 std::cout << "dim(R) = " << dimR << "\n";
156 std::cout << "dim(W) = " << dimW << "\n";
157 std::cout << "dim(R+W) = " << dimR + dimW << "\n";
158 std::cout << "***********************************************************\n";
159 }
160
161 // 7. Define the two BlockStructure of the problem. block_offsets is used
162 // for Vector based on dof (like ParGridFunction or ParLinearForm),
163 // block_trueOffstes is used for Vector based on trueDof (HypreParVector
164 // for the rhs and solution of the linear system). The offsets computed
165 // here are local to the processor.
166 Array<int> block_offsets(3); // number of variables + 1
167 block_offsets[0] = 0;
168 block_offsets[1] = R_space->GetVSize();
169 block_offsets[2] = W_space->GetVSize();
170 block_offsets.PartialSum();
171
172 Array<int> block_trueOffsets(3); // number of variables + 1
173 block_trueOffsets[0] = 0;
174 block_trueOffsets[1] = R_space->TrueVSize();
175 block_trueOffsets[2] = W_space->TrueVSize();
176 block_trueOffsets.PartialSum();
177
178 // 8. Define the coefficients, analytical solution, and rhs of the PDE.
179 ConstantCoefficient k(1.0);
180
184
187
188 // 9. Define the parallel grid function and parallel linear forms, solution
189 // vector and rhs.
190 BlockVector x(block_offsets), rhs(block_offsets);
191 BlockVector trueX(block_trueOffsets), trueRhs(block_trueOffsets);
192
193 ParLinearForm *fform(new ParLinearForm);
194 fform->Update(R_space, rhs.GetBlock(0), 0);
197 fform->Assemble();
198 fform->ParallelAssemble(trueRhs.GetBlock(0));
199
200 ParLinearForm *gform(new ParLinearForm);
201 gform->Update(W_space, rhs.GetBlock(1), 0);
203 gform->Assemble();
204 gform->ParallelAssemble(trueRhs.GetBlock(1));
205
206 // 10. Assemble the finite element matrices for the Darcy operator
207 //
208 // D = [ M B^T ]
209 // [ B 0 ]
210 // where:
211 //
212 // M = \int_\Omega k u_h \cdot v_h d\Omega u_h, v_h \in R_h
213 // B = -\int_\Omega \div u_h q_h d\Omega u_h \in R_h, q_h \in W_h
214 ParBilinearForm *mVarf(new ParBilinearForm(R_space));
215 ParMixedBilinearForm *bVarf(new ParMixedBilinearForm(R_space, W_space));
216
217 PetscParMatrix *pM = NULL, *pB = NULL, *pBT = NULL;
218 HypreParMatrix *M = NULL, *B = NULL, *BT = NULL;
219 Operator::Type tid =
220 !use_petsc ? Operator::Hypre_ParCSR :
221 (use_nonoverlapping ? Operator::PETSC_MATIS : Operator::PETSC_MATAIJ);
222 OperatorHandle Mh(tid), Bh(tid);
223
225 mVarf->Assemble();
226 mVarf->Finalize();
227 mVarf->ParallelAssemble(Mh);
228 if (!use_petsc) { Mh.Get(M); }
229 else { Mh.Get(pM); }
230 Mh.SetOperatorOwner(false);
231
233 bVarf->Assemble();
234 bVarf->Finalize();
235 if (!use_petsc)
236 {
237 B = bVarf->ParallelAssemble();
238 (*B) *= -1;
239 }
240 else
241 {
242 bVarf->ParallelAssemble(Bh);
243 Bh.Get(pB);
244 Bh.SetOperatorOwner(false);
245 (*pB) *= -1;
246 }
247
248 if (!use_petsc) { BT = B->Transpose(); }
249 else { pBT = pB->Transpose(); };
250
251 Operator *darcyOp = NULL;
252 if (!use_petsc)
253 {
254 BlockOperator *tdarcyOp = new BlockOperator(block_trueOffsets);
255 tdarcyOp->SetBlock(0,0,M);
256 tdarcyOp->SetBlock(0,1,BT);
257 tdarcyOp->SetBlock(1,0,B);
258 darcyOp = tdarcyOp;
259 }
260 else
261 {
262 // We construct the BlockOperator and we then convert it to a
263 // PetscParMatrix to avoid any conversion in the construction of the
264 // preconditioners.
265 BlockOperator *tdarcyOp = new BlockOperator(block_trueOffsets);
266 tdarcyOp->SetBlock(0,0,pM);
267 tdarcyOp->SetBlock(0,1,pBT);
268 tdarcyOp->SetBlock(1,0,pB);
269 darcyOp = new PetscParMatrix(pM->GetComm(),tdarcyOp,
270 use_nonoverlapping ? Operator::PETSC_MATIS :
272 delete tdarcyOp;
273 }
274
275 // 11. Construct the operators for preconditioner
276 //
277 // P = [ diag(M) 0 ]
278 // [ 0 B diag(M)^-1 B^T ]
279 //
280 // Here we use Symmetric Gauss-Seidel to approximate the inverse of the
281 // pressure Schur Complement.
282 PetscPreconditioner *pdarcyPr = NULL;
283 BlockDiagonalPreconditioner *darcyPr = NULL;
284 HypreSolver *invM = NULL, *invS = NULL;
285 HypreParMatrix *S = NULL;
286 HypreParMatrix *MinvBt = NULL;
287 HypreParVector *Md = NULL;
288 if (!use_petsc)
289 {
290 MinvBt = B->Transpose();
291 Md = new HypreParVector(MPI_COMM_WORLD, M->GetGlobalNumRows(),
292 M->GetRowStarts());
293 M->GetDiag(*Md);
294
295 MinvBt->InvScaleRows(*Md);
296 S = ParMult(B, MinvBt);
297
298 invM = new HypreDiagScale(*M);
299 invS = new HypreBoomerAMG(*S);
300
301 invM->iterative_mode = false;
302 invS->iterative_mode = false;
303
304 darcyPr = new BlockDiagonalPreconditioner(block_trueOffsets);
305 darcyPr->SetDiagonalBlock(0, invM);
306 darcyPr->SetDiagonalBlock(1, invS);
307 }
308 else
309 {
310 if (use_nonoverlapping)
311 {
313
314 // For saddle point problems, we need to provide BDDC the list of
315 // boundary dofs either essential or natural.
316 // Since R_space is the only space that may have boundary dofs and it
317 // is ordered first then W_space, we don't need any local offset when
318 // specifying the dofs.
319 Array<int> bdr_tdof_list;
320 if (pmesh->bdr_attributes.Size())
321 {
322 Array<int> bdr(pmesh->bdr_attributes.Max());
323 bdr = 1;
324
325 if (!local_bdr_spec)
326 {
327 // Essential dofs in global ordering
328 R_space->GetEssentialTrueDofs(bdr, bdr_tdof_list);
329 }
330 else
331 {
332 // Alternatively, you can also provide the list of dofs in local
333 // ordering
334 R_space->GetEssentialVDofs(bdr, bdr_tdof_list);
335 bdr_tdof_list.SetSize(R_space->GetVSize()+W_space->GetVSize(),0);
336 }
337 opts.SetNatBdrDofs(&bdr_tdof_list,local_bdr_spec);
338 }
339 else
340 {
341 MFEM_WARNING("Missing boundary dofs. This may cause solver failures.");
342 }
343
345 pdarcyPr = new PetscBDDCSolver(MPI_COMM_WORLD,*darcyOp,opts,"prec_");
346 }
347 else
348 {
349 // With PETSc, we can construct the (same) block-diagonal solver with
350 // command line options (see rc_ex5p_fieldsplit)
351 pdarcyPr = new PetscFieldSplitSolver(MPI_COMM_WORLD,*darcyOp,"prec_");
352 }
353 }
354
355 // 12. Solve the linear system with MINRES.
356 // Check the norm of the unpreconditioned residual.
357
358 int maxIter(500);
359 real_t rtol(1.e-6);
360 real_t atol(1.e-10);
361
362 chrono.Clear();
363 chrono.Start();
364
365 trueX = 0.0;
366 if (!use_petsc)
367 {
368 MINRESSolver solver(MPI_COMM_WORLD);
369 solver.SetAbsTol(atol);
370 solver.SetRelTol(rtol);
371 solver.SetMaxIter(maxIter);
372 solver.SetOperator(*darcyOp);
373 solver.SetPreconditioner(*darcyPr);
374 solver.SetPrintLevel(1);
375 solver.Mult(trueRhs, trueX);
376 if (verbose)
377 {
378 if (solver.GetConverged())
379 {
380 std::cout << "MINRES converged in " << solver.GetNumIterations()
381 << " iterations with a residual norm of "
382 << solver.GetFinalNorm() << ".\n";
383 }
384 else
385 {
386 std::cout << "MINRES did not converge in "
387 << solver.GetNumIterations()
388 << " iterations. Residual norm is "
389 << solver.GetFinalNorm() << ".\n";
390 }
391 std::cout << "MINRES solver took " << chrono.RealTime() << "s. \n";
392 }
393
394 }
395 else
396 {
397 std::string solvertype;
398 PetscLinearSolver *solver;
399 if (use_nonoverlapping)
400 {
401 // We can use conjugate gradients to solve the problem
402 solver = new PetscPCGSolver(MPI_COMM_WORLD);
403 solvertype = "PCG";
404 }
405 else
406 {
407 solver = new PetscLinearSolver(MPI_COMM_WORLD);
408 solvertype = "MINRES";
409 }
410 solver->SetOperator(*darcyOp);
411 solver->SetPreconditioner(*pdarcyPr);
412 solver->SetAbsTol(atol);
413 solver->SetRelTol(rtol);
414 solver->SetMaxIter(maxIter);
415 solver->SetPrintLevel(2);
416 solver->Mult(trueRhs, trueX);
417 if (verbose)
418 {
419 if (solver->GetConverged())
420 {
421 std::cout << solvertype << " converged in "
422 << solver->GetNumIterations()
423 << " iterations with a residual norm of "
424 << solver->GetFinalNorm() << ".\n";
425 }
426 else
427 {
428 std::cout << solvertype << " did not converge in "
429 << solver->GetNumIterations()
430 << " iterations. Residual norm is "
431 << solver->GetFinalNorm() << ".\n";
432 }
433 std::cout << solvertype << " solver took "
434 << chrono.RealTime() << "s. \n";
435 }
436 delete solver;
437 }
438 chrono.Stop();
439
440 // 13. Extract the parallel grid function corresponding to the finite element
441 // approximation X. This is the local solution on each processor. Compute
442 // L2 error norms.
445 u->MakeRef(R_space, x.GetBlock(0), 0);
446 p->MakeRef(W_space, x.GetBlock(1), 0);
447 u->Distribute(&(trueX.GetBlock(0)));
448 p->Distribute(&(trueX.GetBlock(1)));
449
450 int order_quad = max(2, 2*order+1);
452 for (int i=0; i < Geometry::NumGeom; ++i)
453 {
455 }
456
457 real_t err_u = u->ComputeL2Error(ucoeff, irs);
458 real_t norm_u = ComputeGlobalLpNorm(2, ucoeff, *pmesh, irs);
459 real_t err_p = p->ComputeL2Error(pcoeff, irs);
460 real_t norm_p = ComputeGlobalLpNorm(2, pcoeff, *pmesh, irs);
461
462 if (verbose)
463 {
464 std::cout << "|| u_h - u_ex || / || u_ex || = " << err_u / norm_u << "\n";
465 std::cout << "|| p_h - p_ex || / || p_ex || = " << err_p / norm_p << "\n";
466 }
467
468 // 14. Save the refined mesh and the solution in parallel. This output can be
469 // viewed later using GLVis: "glvis -np <np> -m mesh -g sol_*".
470 {
471 ostringstream mesh_name, u_name, p_name;
472 mesh_name << "mesh." << setfill('0') << setw(6) << myid;
473 u_name << "sol_u." << setfill('0') << setw(6) << myid;
474 p_name << "sol_p." << setfill('0') << setw(6) << myid;
475
476 ofstream mesh_ofs(mesh_name.str().c_str());
477 mesh_ofs.precision(8);
478 pmesh->Print(mesh_ofs);
479
480 ofstream u_ofs(u_name.str().c_str());
481 u_ofs.precision(8);
482 u->Save(u_ofs);
483
484 ofstream p_ofs(p_name.str().c_str());
485 p_ofs.precision(8);
486 p->Save(p_ofs);
487 }
488
489 // 15. Save data in the VisIt format
490 VisItDataCollection visit_dc("Example5-Parallel", pmesh);
491 visit_dc.RegisterField("velocity", u);
492 visit_dc.RegisterField("pressure", p);
493 visit_dc.SetFormat(!par_format ?
496 visit_dc.Save();
497
498 // 16. Send the solution by socket to a GLVis server.
499 if (visualization)
500 {
501 char vishost[] = "localhost";
502 int visport = 19916;
504 u_sock << "parallel " << num_procs << " " << myid << "\n";
505 u_sock.precision(8);
506 u_sock << "solution\n" << *pmesh << *u << "window_title 'Velocity'"
507 << endl;
508 // Make sure all ranks have sent their 'u' solution before initiating
509 // another set of GLVis connections (one from each rank):
510 MPI_Barrier(pmesh->GetComm());
512 p_sock << "parallel " << num_procs << " " << myid << "\n";
513 p_sock.precision(8);
514 p_sock << "solution\n" << *pmesh << *p << "window_title 'Pressure'"
515 << endl;
516 }
517
518 // 17. Free the used memory.
519 delete fform;
520 delete gform;
521 delete u;
522 delete p;
523 delete darcyOp;
524 delete darcyPr;
525 delete pdarcyPr;
526 delete invM;
527 delete invS;
528 delete S;
529 delete BT;
530 delete B;
531 delete M;
532 delete pBT;
533 delete pB;
534 delete pM;
535 delete MinvBt;
536 delete Md;
537 delete mVarf;
538 delete bVarf;
539 delete W_space;
540 delete R_space;
541 delete l2_coll;
542 delete hdiv_coll;
543 delete pmesh;
544
545 // We finalize PETSc
546 if (use_petsc) { MFEMFinalizePetsc(); }
547
548 return 0;
549}
550
551
552void uFun_ex(const Vector & x, Vector & u)
553{
554 real_t xi(x(0));
555 real_t yi(x(1));
556 real_t zi(0.0);
557 if (x.Size() == 3)
558 {
559 zi = x(2);
560 }
561
562 u(0) = - exp(xi)*sin(yi)*cos(zi);
563 u(1) = - exp(xi)*cos(yi)*cos(zi);
564
565 if (x.Size() == 3)
566 {
567 u(2) = exp(xi)*sin(yi)*sin(zi);
568 }
569}
570
571// Change if needed
573{
574 real_t xi(x(0));
575 real_t yi(x(1));
576 real_t zi(0.0);
577
578 if (x.Size() == 3)
579 {
580 zi = x(2);
581 }
582
583 return exp(xi)*sin(yi)*cos(zi);
584}
585
586void fFun(const Vector & x, Vector & f)
587{
588 f = 0.0;
589}
590
592{
593 if (x.Size() == 3)
594 {
595 return -pFun_ex(x);
596 }
597 else
598 {
599 return 0;
600 }
601}
602
604{
605 return (-pFun_ex(x));
606}
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition array.hpp:697
int Size() const
Return the logical size of the array.
Definition array.hpp:144
void PartialSum()
Fill the entries of the array with the cumulative sum of the entries.
Definition array.cpp:103
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 SetFormat(int fmt)
Set the desired output mesh and data format.
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
Abstract class for hypre's solvers and preconditioners.
Definition hypre.hpp:1162
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
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition mesh.hpp:282
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)
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
void SetOperatorOwner(bool own=true)
Set the ownership flag for the held Operator.
Definition handle.hpp:120
void Get(OpType *&A) const
Return the Operator pointer statically cast to a given OpType.
Definition handle.hpp:114
Abstract operator.
Definition operator.hpp:25
Type
Enumeration defining IDs for some classes derived from Operator.
Definition operator.hpp:284
@ PETSC_MATIS
ID for class PetscParMatrix, MATIS format.
Definition operator.hpp:289
@ Hypre_ParCSR
ID for class HypreParMatrix.
Definition operator.hpp:287
@ PETSC_MATAIJ
ID for class PetscParMatrix, MATAIJ format.
Definition operator.hpp:288
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.
Abstract parallel finite element space.
Definition pfespace.hpp:29
void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1) const override
HYPRE_BigInt GlobalTrueVSize() const
Definition pfespace.hpp:285
void GetEssentialVDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_dofs, int component=-1) const override
Determine the boundary degrees of freedom.
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.
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P_test^t A P_trial.
Auxiliary class for BDDC customization.
Definition petsc.hpp:831
void SetNatBdrDofs(const Array< int > *natdofs, bool loc=false)
Specify dofs on the natural boundary.
Definition petsc.hpp:858
Abstract class for PETSc's linear solvers.
Definition petsc.hpp:750
virtual void SetOperator(const Operator &op)
Definition petsc.cpp:2965
void SetPreconditioner(Solver &precond)
Definition petsc.cpp:3110
virtual void Mult(const Vector &b, Vector &x) const
Application of the solver.
Definition petsc.cpp:3190
Wrapper for PETSc's matrix class.
Definition petsc.hpp:319
MPI_Comm GetComm() const
Get the associated MPI communicator.
Definition petsc.cpp:1330
Abstract class for PETSc's preconditioners.
Definition petsc.hpp:805
void SetAbsTol(real_t tol)
Definition petsc.cpp:2401
void SetPrintLevel(int plev)
Definition petsc.cpp:2453
void SetMaxIter(int max_iter)
Definition petsc.cpp:2426
void SetRelTol(real_t tol)
Definition petsc.cpp:2376
int GetNumIterations()
Definition petsc.cpp:2694
real_t GetFinalNorm()
Definition petsc.cpp:2727
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition fe_coll.hpp:386
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
for VectorFiniteElements (Nedelec, Raviart-Thomas)
Definition lininteg.hpp:347
A general vector function coefficient.
Vector data type.
Definition vector.hpp:80
int Size() const
Returns the size of the vector.
Definition vector.hpp:218
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
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. .
void MFEMInitializePetsc()
Convenience functions to initialize/finalize PETSc.
Definition petsc.cpp:201
void MFEMFinalizePetsc()
Definition petsc.cpp:241
HypreParMatrix * ParMult(const HypreParMatrix *A, const HypreParMatrix *B, bool own_matrix)
Definition hypre.cpp:2949
float real_t
Definition config.hpp:43
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30
A global object with all integration rules (defined in intrules.cpp)
Definition intrules.hpp:486
const char vishost[]
real_t p(const Vector &x, real_t t)
real_t pFun_ex(const Vector &x)
Definition ex5p.cpp:572
void uFun_ex(const Vector &x, Vector &u)
Definition ex5p.cpp:552
real_t gFun(const Vector &x)
Definition ex5p.cpp:591
real_t f_natural(const Vector &x)
Definition ex5p.cpp:603
void fFun(const Vector &x, Vector &f)
Definition ex5p.cpp:586