MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
ex1p.cpp
Go to the documentation of this file.
1 // MFEM Example 1 - Parallel Version
2 // PUMI Modification
3 //
4 // Compile with: make ex1p
5 //
6 // Sample runs:
7 // mpirun -np 8 ex1p -m ../../data/pumi/parallel/Kova/Kova100k_8.smb
8 // -p ../../data/pumi/geom/Kova.dmg -o 1 -go 2
9 //
10 // Note: Example models + meshes for the PUMI examples can be downloaded
11 // from github.com/mfem/data/pumi. After downloading we recommend
12 // creating a symbolic link to the above directory in ../../data.
13 //
14 // Description: This example code demonstrates the use of MFEM to define a
15 // simple finite element discretization of the Laplace problem
16 // -Delta u = 1 with homogeneous Dirichlet boundary conditions.
17 // Specifically, we discretize using a FE space of the specified
18 // order, or if order < 1 using an isoparametric/isogeometric
19 // space (i.e. quadratic for quadratic curvilinear mesh, NURBS for
20 // NURBS mesh, etc.)
21 //
22 // The example highlights the use of mesh refinement, finite
23 // element grid functions, as well as linear and bilinear forms
24 // corresponding to the left-hand side and right-hand side of the
25 // discrete linear system. We also cover the explicit elimination
26 // of essential boundary conditions, static condensation, and the
27 // optional connection to the GLVis tool for visualization.
28 //
29 // This PUMI modification demonstrates how PUMI's API can be used
30 // to load a parallel PUMI mesh classified on a geometric model
31 // and then generate the corresponding parallel MFEM mesh. The
32 // example also performs a "uniform" refinement, similar to the
33 // MFEM examples, for coarse meshes. However, the refinement is
34 // performed using the PUMI API. The inputs are a Parasolid
35 // model, "*.xmt_txt" and SCOREC parallel meshes "*.smb". The
36 // option "-o" is used for the Finite Element order and "-go" for
37 // the geometry order. Note that they can be used independently:
38 // "-o 8 -go 3" solves for 8th order FE on third order geometry.
39 //
40 // NOTE: Model/Mesh files for this example are in the (large) data file
41 // repository of MFEM here https://github.com/mfem/data under the
42 // folder named "pumi", which consists of the following sub-folders:
43 // a) geom --> model files
44 // b) parallel --> parallel pumi mesh files
45 // c) serial --> serial pumi mesh files
46 
47 
48 #include "mfem.hpp"
49 #include <fstream>
50 #include <iostream>
51 
52 #ifdef MFEM_USE_SIMMETRIX
53 #include <SimUtil.h>
54 #include <gmi_sim.h>
55 #endif
56 #include <apfMDS.h>
57 #include <gmi_null.h>
58 #include <PCU.h>
59 #include <apfConvert.h>
60 #include <gmi_mesh.h>
61 #include <crv.h>
62 
63 #ifndef MFEM_USE_PUMI
64 #error This example requires that MFEM is built with MFEM_USE_PUMI=YES
65 #endif
66 
67 using namespace std;
68 using namespace mfem;
69 
70 int main(int argc, char *argv[])
71 {
72  // 1. Initialize MPI and HYPRE.
73  Mpi::Init(argc, argv);
74  int num_procs = Mpi::WorldSize();
75  int myid = Mpi::WorldRank();
76  Hypre::Init();
77 
78  // 2. Parse command-line options.
79  const char *mesh_file = "../../data/pumi/parallel/Kova/Kova100k_8.smb";
80 #ifdef MFEM_USE_SIMMETRIX
81  const char *model_file = "../../data/pumi/geom/Kova.x_t";
82 #else
83  const char *model_file = "../../data/pumi/geom/Kova.dmg";
84 #endif
85  int order = 1;
86  bool static_cond = false;
87  bool visualization = 1;
88  int geom_order = 1;
89 
90  OptionsParser args(argc, argv);
91  args.AddOption(&mesh_file, "-m", "--mesh",
92  "Mesh file to use.");
93  args.AddOption(&order, "-o", "--order",
94  "Finite element order (polynomial degree) or -1 for"
95  " isoparametric space.");
96  args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
97  "--no-static-condensation", "Enable static condensation.");
98  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
99  "--no-visualization",
100  "Enable or disable GLVis visualization.");
101  args.AddOption(&model_file, "-p", "--parasolid",
102  "Parasolid model to use.");
103  args.AddOption(&geom_order, "-go", "--geometry_order",
104  "Geometric order of the model");
105  args.Parse();
106  if (!args.Good())
107  {
108  if (myid == 0)
109  {
110  args.PrintUsage(cout);
111  }
112  return 1;
113  }
114  if (myid == 0)
115  {
116  args.PrintOptions(cout);
117  }
118 
119  // 3. Read the SCOREC Mesh
120  PCU_Comm_Init();
121 #ifdef MFEM_USE_SIMMETRIX
122  Sim_readLicenseFile(0);
123  gmi_sim_start();
124  gmi_register_sim();
125 #endif
126  gmi_register_mesh();
127 
128  apf::Mesh2* pumi_mesh;
129  pumi_mesh = apf::loadMdsMesh(model_file, mesh_file);
130 
131  // 4. Increase the geometry order and refine the mesh if necessary. Parallel
132  // uniform refinement is performed if the total number of elements is less
133  // than 10,000.
134  int dim = pumi_mesh->getDimension();
135  int nEle = pumi_mesh->count(dim);
136  int ref_levels = (int)floor(log(10000./nEle)/log(2.)/dim);
137 
138  if (geom_order > 1)
139  {
140  crv::BezierCurver bc(pumi_mesh, geom_order, 2);
141  bc.run();
142  }
143 
144  // Perform Uniform refinement
145  if (ref_levels > 1)
146  {
147  auto uniInput = ma::configureUniformRefine(pumi_mesh, ref_levels);
148 
149  if (geom_order > 1)
150  {
151  crv::adapt(uniInput);
152  }
153  else
154  {
155  ma::adapt(uniInput);
156  }
157  }
158 
159  pumi_mesh->verify();
160 
161  // 5. Create the parallel MFEM mesh object from the parallel PUMI mesh.
162  // We can handle triangular and tetrahedral meshes. Note that the
163  // mesh resolution is performed on the PUMI mesh.
164  ParMesh *pmesh = new ParPumiMesh(MPI_COMM_WORLD, pumi_mesh);
165 
166  // 6. Define a parallel finite element space on the parallel mesh. Here we
167  // use continuous Lagrange finite elements of the specified order. If
168  // order < 1, we instead use an isoparametric/isogeometric space.
170  if (order > 0)
171  {
172  fec = new H1_FECollection(order, dim);
173  }
174  else if (pmesh->GetNodes())
175  {
176  fec = pmesh->GetNodes()->OwnFEC();
177  if (myid == 0)
178  {
179  cout << "Using isoparametric FEs: " << fec->Name() << endl;
180  }
181  }
182  else
183  {
184  fec = new H1_FECollection(order = 1, dim);
185  }
186  ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
187  HYPRE_BigInt size = fespace->GlobalTrueVSize();
188  if (myid == 0)
189  {
190  cout << "Number of finite element unknowns: " << size << endl;
191  }
192 
193  // 7. Determine the list of true (i.e. parallel conforming) essential
194  // boundary dofs. In this example, the boundary conditions are defined
195  // by marking all the boundary attributes from the mesh as essential
196  // (Dirichlet) and converting them to a list of true dofs.
197  Array<int> ess_tdof_list;
198  if (pmesh->bdr_attributes.Size())
199  {
200  Array<int> ess_bdr(pmesh->bdr_attributes.Max());
201  ess_bdr = 1;
202  fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
203  }
204 
205  // 8. Set up the parallel linear form b(.) which corresponds to the
206  // right-hand side of the FEM linear system, which in this case is
207  // (1,phi_i) where phi_i are the basis functions in fespace.
208  ParLinearForm *b = new ParLinearForm(fespace);
209  ConstantCoefficient one(1.0);
211  b->Assemble();
212 
213  // 9. Define the solution vector x as a parallel finite element grid function
214  // corresponding to fespace. Initialize x with initial guess of zero,
215  // which satisfies the boundary conditions.
216  ParGridFunction x(fespace);
217  x = 0.0;
218 
219  // 10. Set up the parallel bilinear form a(.,.) on the finite element space
220  // corresponding to the Laplacian operator -Delta, by adding the Diffusion
221  // domain integrator.
222  ParBilinearForm *a = new ParBilinearForm(fespace);
224 
225  // 11. Assemble the parallel bilinear form and the corresponding linear
226  // system, applying any necessary transformations such as: parallel
227  // assembly, eliminating boundary conditions, applying conforming
228  // constraints for non-conforming AMR, static condensation, etc.
229  if (static_cond) { a->EnableStaticCondensation(); }
230  a->Assemble();
231 
232  HypreParMatrix A;
233  Vector B, X;
234  a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
235 
236  if (myid == 0)
237  {
238  cout << "Size of linear system: " << A.GetGlobalNumRows() << endl;
239  }
240 
241  // 12. Define and apply a parallel PCG solver for AX=B with the BoomerAMG
242  // preconditioner from hypre.
243  HypreSolver *amg = new HypreBoomerAMG(A);
244  HyprePCG *pcg = new HyprePCG(A);
245  pcg->SetTol(1e-12);
246  pcg->SetMaxIter(200);
247  pcg->SetPrintLevel(2);
248  pcg->SetPreconditioner(*amg);
249  pcg->Mult(B, X);
250 
251  // 13. Recover the parallel grid function corresponding to X. This is the
252  // local finite element solution on each processor.
253  a->RecoverFEMSolution(X, *b, x);
254 
255  // 14. Save the refined mesh and the solution in parallel. This output can
256  // be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
257  {
258  ostringstream mesh_name, sol_name;
259  mesh_name << "mesh." << setfill('0') << setw(6) << myid;
260  sol_name << "sol." << setfill('0') << setw(6) << myid;
261 
262  ofstream mesh_ofs(mesh_name.str().c_str());
263  mesh_ofs.precision(8);
264  pmesh->Print(mesh_ofs);
265 
266  ofstream sol_ofs(sol_name.str().c_str());
267  sol_ofs.precision(8);
268  x.Save(sol_ofs);
269  }
270 
271  // 15. Send the solution by socket to a GLVis server.
272  if (visualization)
273  {
274  char vishost[] = "localhost";
275  int visport = 19916;
276  socketstream sol_sock(vishost, visport);
277  sol_sock << "parallel " << num_procs << " " << myid << "\n";
278  sol_sock.precision(8);
279  sol_sock << "solution\n" << *pmesh << x << flush;
280  }
281 
282  // 16. Free the used memory.
283  delete pcg;
284  delete amg;
285  delete a;
286  delete b;
287  delete fespace;
288  if (order > 0) { delete fec; }
289  delete pmesh;
290 
291  pumi_mesh->destroyNative();
292  apf::destroyMesh(pumi_mesh);
293  PCU_Comm_Free();
294 
295 #ifdef MFEM_USE_SIMMETRIX
296  gmi_sim_stop();
297  Sim_unregisterAllKeys();
298 #endif
299 
300  return 0;
301 }
void SetTol(double tol)
Definition: hypre.cpp:3995
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:108
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Definition: pfespace.cpp:1032
A coefficient that is constant across space and time.
Definition: coefficient.hpp:84
virtual void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B, int copy_interior=0)
Form the linear system A X = B, corresponding to this bilinear form and the linear form b(...
HYPRE_BigInt GlobalTrueVSize() const
Definition: pfespace.hpp:285
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:873
Abstract parallel finite element space.
Definition: pfespace.hpp:28
void SetPrintLevel(int print_lvl)
Definition: hypre.cpp:4015
Class for PUMI parallel meshes.
Definition: pumi.hpp:69
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1579
Class for parallel linear form.
Definition: plinearform.hpp:26
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
constexpr char vishost[]
void Assemble()
Assembles the ParLinearForm i.e. sums over all domain/bdr integrators.
Definition: plinearform.cpp:46
double b
Definition: lissajous.cpp:42
constexpr int visport
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)
Assemble the local matrix.
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
void SetMaxIter(int max_iter)
Definition: hypre.cpp:4005
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:41
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:270
PCG solver in hypre.
Definition: hypre.hpp:1204
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
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
double a
Definition: lissajous.cpp:41
virtual const char * Name() const
Definition: fe_coll.hpp:65
HYPRE_BigInt GetGlobalNumRows() const
Return the global number of rows.
Definition: hypre.hpp:635
int dim
Definition: ex24.cpp:53
void SetPreconditioner(HypreSolver &precond)
Set the hypre solver to be used as a preconditioner.
Definition: hypre.cpp:4020
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
Class for parallel bilinear form.
Abstract class for hypre&#39;s solvers and preconditioners.
Definition: hypre.hpp:1092
Vector data type.
Definition: vector.hpp:60
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:7908
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:220
void Print(std::ostream &out=mfem::out) const override
Definition: pmesh.cpp:4770
Class for parallel grid function.
Definition: pgridfunc.hpp:32
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Solve Ax=b with hypre&#39;s PCG.
Definition: hypre.cpp:4043
Class for parallel meshes.
Definition: pmesh.hpp:32
int main()
void EnableStaticCondensation()
Enable the use of static condensation. For details see the description for class StaticCondensation i...
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150