MFEM  v4.5.0
Finite element discretization library
ex6p.cpp
Go to the documentation of this file.
1 // MFEM Example 6 - Parallel Version
2 // PUMI Modification
3 //
4 // Compile with: make ex6p
5 //
6 // Sample runs: mpirun -np 8 ex6p
7 //
8 // Description: This is a version of Example 1 with a simple adaptive mesh
9 // refinement loop. The problem being solved is again the Laplace
10 // equation -Delta u = 1 with homogeneous Dirichlet boundary
11 // conditions. The problem is solved on a sequence of meshes which
12 // are adapted in a conforming (tetrahedrons) manner according
13 // to a simple SPR ZZ error estimator.
14 //
15 // This PUMI variation also performs a "uniform" refinement,
16 // similar to MFEM examples, for coarse meshes. However, the
17 // refinement is performed using the PUMI API. A new option "-ar"
18 // is added to modify the "adapt_ratio" which is the fraction of
19 // allowable error that scales the output size field of the error
20 // estimator.
21 //
22 // NOTE: Model/Mesh files for this example are in the (large) data file
23 // repository of MFEM here https://github.com/mfem/data under the
24 // folder named "pumi", which consists of the following sub-folders:
25 // a) geom --> model files
26 // b) parallel --> parallel pumi mesh files
27 // c) serial --> serial pumi mesh files
28 
29 #include "mfem.hpp"
30 #include <fstream>
31 #include <iostream>
32 
33 #ifdef MFEM_USE_SIMMETRIX
34 #include <SimUtil.h>
35 #include <gmi_sim.h>
36 #endif
37 #include <apfMDS.h>
38 #include <gmi_null.h>
39 #include <PCU.h>
40 #include <spr.h>
41 #include <apfConvert.h>
42 #include <gmi_mesh.h>
43 #include <crv.h>
44 
45 #ifndef MFEM_USE_PUMI
46 #error This example requires that MFEM is built with MFEM_USE_PUMI=YES
47 #endif
48 
49 using namespace std;
50 using namespace mfem;
51 
52 int main(int argc, char *argv[])
53 {
54  // 1. Initialize MPI and HYPRE.
55  Mpi::Init(argc, argv);
56  int num_procs = Mpi::WorldSize();
57  int myid = Mpi::WorldRank();
58  Hypre::Init();
59 
60  // 2. Parse command-line options.
61  const char *mesh_file = "../../data/pumi/parallel/Kova/Kova100k_8.smb";
62 #ifdef MFEM_USE_SIMMETRIX
63  const char *model_file = "../../data/pumi/geom/Kova.x_t";
64  const char *smd_file = NULL;
65 #else
66  const char *model_file = "../../data/pumi/geom/Kova.dmg";
67 #endif
68  int order = 1;
69  bool static_cond = false;
70  bool visualization = 1;
71  int geom_order = 1;
72  double adapt_ratio = 0.05;
73 
74  OptionsParser args(argc, argv);
75  args.AddOption(&mesh_file, "-m", "--mesh",
76  "Mesh file to use.");
77  args.AddOption(&order, "-o", "--order",
78  "Finite element order (polynomial degree) or -1 for"
79  " isoparametric space.");
80  args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
81  "--no-static-condensation", "Enable static condensation.");
82  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
83  "--no-visualization",
84  "Enable or disable GLVis visualization.");
85  args.AddOption(&model_file, "-p", "--model",
86  "parasolid or .dmg model to use.");
87 #ifdef MFEM_USE_SIMMETRIX
88  args.AddOption(&smd_file, "-sm", "--smd_model",
89  "smd model file to use.");
90 #endif
91  args.AddOption(&geom_order, "-go", "--geometry_order",
92  "Geometric order of the model");
93  args.AddOption(&adapt_ratio, "-ar", "--adapt_ratio",
94  "adaptation factor used in MeshAdapt");
95  args.Parse();
96  if (!args.Good())
97  {
98  if (myid == 0)
99  {
100  args.PrintUsage(cout);
101  }
102  return 1;
103  }
104  if (myid == 0)
105  {
106  args.PrintOptions(cout);
107  }
108 
109  // 3. Read the SCOREC Mesh.
110  PCU_Comm_Init();
111 #ifdef MFEM_USE_SIMMETRIX
112  Sim_readLicenseFile(0);
113  gmi_sim_start();
114  gmi_register_sim();
115 #endif
116  gmi_register_mesh();
117 
118  apf::Mesh2* pumi_mesh;
119 #ifdef MFEM_USE_SIMMETRIX
120  if (smd_file)
121  {
122  gmi_model *mixed_model = gmi_sim_load(model_file, smd_file);
123  pumi_mesh = apf::loadMdsMesh(mixed_model, mesh_file);
124  }
125  else
126 #endif
127  {
128  pumi_mesh = apf::loadMdsMesh(model_file, mesh_file);
129  }
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 100,000.
134  int dim = pumi_mesh->getDimension();
135  int nEle = pumi_mesh->count(dim);
136  int ref_levels = (int)floor(log(100000./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 (myid == 1)
146  {
147  std::cout << " ref level : " << ref_levels << std::endl;
148  }
149 
150  if (ref_levels > 1)
151  {
152  auto uniInput = ma::configureUniformRefine(pumi_mesh, ref_levels);
153 
154  if ( geom_order > 1)
155  {
156  crv::adapt(uniInput);
157  }
158  else
159  {
160  ma::adapt(uniInput);
161  }
162  }
163 
164  pumi_mesh->verify();
165 
166  // 5. Create the parallel MFEM mesh object from the parallel PUMI mesh. We
167  // can handle triangular and tetrahedral meshes. Note that the mesh
168  // resolution is performed on the PUMI mesh.
169  ParMesh *pmesh = new ParPumiMesh(MPI_COMM_WORLD, pumi_mesh);
170 
171  // 6. Define a parallel finite element space on the parallel mesh. Here we
172  // use continuous Lagrange finite elements of the specified order. If
173  // order < 1, we instead use an isoparametric/isogeometric space.
175  if (order > 0)
176  {
177  fec = new H1_FECollection(order, dim);
178  }
179  else if (pmesh->GetNodes())
180  {
181  fec = pmesh->GetNodes()->OwnFEC();
182  if (myid == 1)
183  {
184  cout << "Using isoparametric FEs: " << fec->Name() << endl;
185  }
186  }
187  else
188  {
189  fec = new H1_FECollection(order = 1, dim);
190  }
191  ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
192  HYPRE_BigInt size = fespace->GlobalTrueVSize();
193  if (myid == 1)
194  {
195  cout << "Number of finite element unknowns: " << size << endl;
196  }
197 
198  // 7. Set up the parallel linear form b(.) which corresponds to the
199  // right-hand side of the FEM linear system, which in this case is
200  // (1,phi_i) where phi_i are the basis functions in fespace.
201  ParLinearForm *b = new ParLinearForm(fespace);
202  ConstantCoefficient one(1.0);
203  b->AddDomainIntegrator(new DomainLFIntegrator(one));
204 
205  // 8. Define the solution vector x as a parallel finite element grid function
206  // corresponding to fespace. Initialize x with initial guess of zero,
207  // which satisfies the boundary conditions.
208  ParGridFunction x(fespace);
209  x = 0.0;
210 
211  // 9. Connect to GLVis.
212  char vishost[] = "localhost";
213  int visport = 19916;
214 
215  socketstream sout;
216  if (visualization)
217  {
218  sout.open(vishost, visport);
219  if (!sout)
220  {
221  if (myid == 0)
222  {
223  cout << "Unable to connect to GLVis server at "
224  << vishost << ':' << visport << endl;
225  cout << "GLVis visualization disabled.\n";
226  }
227  visualization = false;
228  }
229 
230  sout.precision(8);
231  }
232 
233  // 10. Set up the parallel bilinear form a(.,.) on the finite element space
234  // corresponding to the Laplacian operator -Delta, by adding the
235  // Diffusion domain integrator.
236  ParBilinearForm *a = new ParBilinearForm(fespace);
237  a->AddDomainIntegrator(new DiffusionIntegrator(one));
238 
239  // 11. Assemble the parallel bilinear form and the corresponding linear
240  // system, applying any necessary transformations such as: parallel
241  // assembly, eliminating boundary conditions, applying conforming
242  // constraints for non-conforming AMR, static condensation, etc.
243  if (static_cond) { a->EnableStaticCondensation(); }
244 
245  // 12. The main AMR loop. In each iteration we solve the problem on the
246  // current mesh, visualize the solution, and adapt the mesh.
247  apf::Field* Tmag_field = 0;
248  apf::Field* temp_field = 0;
249  apf::Field* ipfield = 0;
250  apf::Field* sizefield = 0;
251  int max_iter = 3;
252 
253  for (int Itr = 0; Itr < max_iter; Itr++)
254  {
255  HYPRE_BigInt global_dofs = fespace->GlobalTrueVSize();
256  if (myid == 1)
257  {
258  cout << "\nAMR iteration " << Itr << endl;
259  cout << "Number of unknowns: " << global_dofs << endl;
260  }
261 
262  // Assemble.
263  a->Assemble();
264  b->Assemble();
265 
266  // Essential boundary condition.
267  Array<int> ess_tdof_list;
268  if (pmesh->bdr_attributes.Size())
269  {
270  Array<int> ess_bdr(pmesh->bdr_attributes.Max());
271  ess_bdr = 1;
272  fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
273  }
274 
275  // Form linear system.
276  HypreParMatrix A;
277  Vector B, X;
278  const int copy_interior = 1;
279  a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B, copy_interior);
280 
281  // 13. Define and apply a parallel PCG solver for AX=B with the BoomerAMG
282  // preconditioner from hypre.
283  HypreBoomerAMG amg;
284  amg.SetPrintLevel(0);
285  CGSolver pcg(A.GetComm());
286  pcg.SetPreconditioner(amg);
287  pcg.SetOperator(A);
288  pcg.SetRelTol(1e-6);
289  pcg.SetMaxIter(200);
290  pcg.SetPrintLevel(3); // print the first and the last iterations only
291  pcg.Mult(B, X);
292 
293  // 14. Recover the parallel grid function corresponding to X. This is the
294  // local finite element solution on each processor.
295  a->RecoverFEMSolution(X, *b, x);
296 
297  // 15. Save in parallel the displaced mesh and the inverted solution (which
298  // gives the backward displacements to the original grid). This output
299  // can be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
300  {
301  ostringstream mesh_name, sol_name;
302  mesh_name << "mesh." << setfill('0') << setw(6) << myid;
303  sol_name << "sol." << setfill('0') << setw(6) << myid;
304 
305  ofstream mesh_ofs(mesh_name.str().c_str());
306  mesh_ofs.precision(8);
307  pmesh->Print(mesh_ofs);
308 
309  ofstream sol_ofs(sol_name.str().c_str());
310  sol_ofs.precision(8);
311  x.Save(sol_ofs);
312  }
313 
314  // 16. Send the above data by socket to a GLVis server. Use the "n" and "b"
315  // keys in GLVis to visualize the displacements.
316  if (visualization)
317  {
318  sout << "parallel " << num_procs << " " << myid << "\n";
319  sout << "solution\n" << *pmesh << x << flush;
320  }
321 
322  // 17. Field transfer. Scalar solution field and magnitude field for error
323  // estimation are created the PUMI mesh.
324  if (order > geom_order)
325  {
326  Tmag_field = apf::createField(pumi_mesh, "field_mag",
327  apf::SCALAR, apf::getLagrange(order));
328  temp_field = apf::createField(pumi_mesh, "T_field",
329  apf::SCALAR, apf::getLagrange(order));
330  }
331  else
332  {
333  Tmag_field = apf::createFieldOn(pumi_mesh, "field_mag",apf::SCALAR);
334  temp_field = apf::createFieldOn(pumi_mesh, "T_field", apf::SCALAR);
335  }
336 
337  ParPumiMesh* pPPmesh = dynamic_cast<ParPumiMesh*>(pmesh);
338  pPPmesh->FieldMFEMtoPUMI(pumi_mesh, &x, temp_field, Tmag_field);
339 
340  ipfield= spr::getGradIPField(Tmag_field, "MFEM_gradip", 2);
341  sizefield = spr::getSPRSizeField(ipfield, adapt_ratio);
342 
343  apf::destroyField(Tmag_field);
344  apf::destroyField(ipfield);
345 
346  // 18. Perform MesAdapt.
347  auto erinput = ma::configure(pumi_mesh, sizefield);
348  if ( geom_order > 1)
349  {
350  crv::adapt(erinput);
351  }
352  else
353  {
354  ma::adapt(erinput);
355  }
356 
357  ParMesh* Adapmesh = new ParPumiMesh(MPI_COMM_WORLD, pumi_mesh);
358  pPPmesh->UpdateMesh(Adapmesh);
359  delete Adapmesh;
360 
361  // 19. Update the FiniteElementSpace, GridFunction, and bilinear form.
362  fespace->Update();
363  x.Update();
364  x = 0.0;
365 
366  pPPmesh->FieldPUMItoMFEM(pumi_mesh, temp_field, &x);
367  a->Update();
368  b->Update();
369 
370  // Destroy fields.
371  apf::destroyField(temp_field);
372  apf::destroyField(sizefield);
373  }
374 
375  // 20. Free the used memory.
376  delete a;
377  delete b;
378  delete fespace;
379  if (order > 0) { delete fec; }
380  delete pmesh;
381 
382  pumi_mesh->destroyNative();
383  apf::destroyMesh(pumi_mesh);
384  PCU_Comm_Free();
385 
386 #ifdef MFEM_USE_SIMMETRIX
387  gmi_sim_stop();
388  Sim_unregisterAllKeys();
389 #endif
390 
391  return 0;
392 }
void UpdateMesh(const ParMesh *AdaptedpMesh)
Update the mesh after adaptation.
Definition: pumi.cpp:956
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
Conjugate gradient method.
Definition: solvers.hpp:465
A coefficient that is constant across space and time.
Definition: coefficient.hpp:84
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
virtual void Update(bool want_transform=true)
Definition: pfespace.cpp:3287
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:68
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150
Abstract parallel finite element space.
Definition: pfespace.hpp:28
void FieldMFEMtoPUMI(apf::Mesh2 *apf_mesh, ParGridFunction *grid_vel, ParGridFunction *grid_pr, apf::Field *vel_field, apf::Field *pr_field, apf::Field *vel_mag_field)
Transfer field from MFEM mesh to PUMI mesh [Mixed].
Definition: pumi.cpp:1216
STL namespace.
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[]
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Definition: pgridfunc.cpp:81
double b
Definition: lissajous.cpp:42
void SetPrintLevel(int print_level)
Definition: hypre.hpp:1659
constexpr int visport
virtual const char * Name() const
Definition: fe_coll.hpp:65
HYPRE_BigInt GlobalTrueVSize() const
Definition: pfespace.hpp:285
MPI_Comm GetComm() const
MPI communicator.
Definition: hypre.hpp:534
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:270
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
void FieldPUMItoMFEM(apf::Mesh2 *apf_mesh, apf::Field *field, ParGridFunction *grid)
Transfer a field from PUMI to MFEM after mesh adapt [Scalar and Vector].
Definition: pumi.cpp:1417
HYPRE_Int HYPRE_BigInt
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:873
double a
Definition: lissajous.cpp:41
int dim
Definition: ex24.cpp:53
Class for parallel bilinear form.
int open(const char hostname[], int port)
Open the socket stream on &#39;port&#39; at &#39;hostname&#39;.
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
Vector data type.
Definition: vector.hpp:60
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:173
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:220
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:7894
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
Class for parallel meshes.
Definition: pmesh.hpp:32
int main(int argc, char *argv[])
Definition: ex6p.cpp:53