1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
|
// MFEM Example 0 - Parallel Version
//
// Compile with: make ex0p
//
// Sample runs: mpirun -np 4 ex0p
// mpirun -np 4 ex0p -m ../data/fichera.mesh
// mpirun -np 4 ex0p -m ../data/square-disc.mesh -o 2
//
// Description: This example code demonstrates the most basic parallel usage of
// MFEM to define a simple finite element discretization of the
// Laplace problem -Delta u = 1 with zero Dirichlet boundary
// conditions. General 2D/3D serial mesh files and finite element
// polynomial degrees can be specified by command line options.
#include "mfem.hpp"
#include <fstream>
#include <iostream>
using namespace std;
using namespace mfem;
int main(int argc, char *argv[])
{
// 1. Initialize MPI and HYPRE.
Mpi::Init(argc, argv);
Hypre::Init();
// 2. Parse command line options.
string mesh_file = "../data/star.mesh";
int order = 1;
OptionsParser args(argc, argv);
args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
args.AddOption(&order, "-o", "--order", "Finite element polynomial degree");
args.ParseCheck();
// 3. Read the serial mesh from the given mesh file.
Mesh serial_mesh(mesh_file);
// 4. Define a parallel mesh by a partitioning of the serial mesh. Refine
// this mesh once in parallel to increase the resolution.
ParMesh mesh(MPI_COMM_WORLD, serial_mesh);
serial_mesh.Clear(); // the serial mesh is no longer needed
mesh.UniformRefinement();
// 5. Define a finite element space on the mesh. Here we use H1 continuous
// high-order Lagrange finite elements of the given order.
H1_FECollection fec(order, mesh.Dimension());
ParFiniteElementSpace fespace(&mesh, &fec);
HYPRE_BigInt total_num_dofs = fespace.GlobalTrueVSize();
if (Mpi::Root())
{
cout << "Number of unknowns: " << total_num_dofs << endl;
}
// 6. Extract the list of all the boundary DOFs. These will be marked as
// Dirichlet in order to enforce zero boundary conditions.
Array<int> boundary_dofs;
fespace.GetBoundaryTrueDofs(boundary_dofs);
// 7. Define the solution x as a finite element grid function in fespace. Set
// the initial guess to zero, which also sets the boundary conditions.
ParGridFunction x(&fespace);
x = 0.0;
// 8. Set up the linear form b(.) corresponding to the right-hand side.
ConstantCoefficient one(1.0);
ParLinearForm b(&fespace);
b.AddDomainIntegrator(new DomainLFIntegrator(one));
b.Assemble();
// 9. Set up the bilinear form a(.,.) corresponding to the -Delta operator.
ParBilinearForm a(&fespace);
a.AddDomainIntegrator(new DiffusionIntegrator);
a.Assemble();
// 10. Form the linear system A X = B. This includes eliminating boundary
// conditions, applying AMR constraints, parallel assembly, etc.
HypreParMatrix A;
Vector B, X;
a.FormLinearSystem(boundary_dofs, x, b, A, X, B);
// 11. Solve the system using PCG with hypre's BoomerAMG preconditioner.
HypreBoomerAMG M(A);
CGSolver cg(MPI_COMM_WORLD);
cg.SetRelTol(1e-12);
cg.SetMaxIter(2000);
cg.SetPrintLevel(1);
cg.SetPreconditioner(M);
cg.SetOperator(A);
cg.Mult(B, X);
// 12. Recover the solution x as a grid function and save to file. The output
// can be viewed using GLVis as follows: "glvis -np <np> -m mesh -g sol"
a.RecoverFEMSolution(X, b, x);
x.Save("sol");
mesh.Save("mesh");
return 0;
}
|
