MFEM  v4.1.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
mesh-explorer.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 //
12 // -----------------------------------------------------
13 // Mesh Explorer Miniapp: Explore and manipulate meshes
14 // -----------------------------------------------------
15 //
16 // This miniapp is a handy tool to examine, visualize and manipulate a given
17 // mesh. Some of its features are:
18 //
19 // - visualizing of mesh materials and individual mesh elements
20 // - mesh scaling, randomization, and general transformation
21 // - manipulation of the mesh curvature
22 // - the ability to simulate parallel partitioning
23 // - quantitative and visual reports of mesh quality
24 //
25 // Compile with: make mesh-explorer
26 //
27 // Sample runs: mesh-explorer
28 // mesh-explorer -m ../../data/beam-tri.mesh
29 // mesh-explorer -m ../../data/star-q2.mesh
30 // mesh-explorer -m ../../data/disc-nurbs.mesh
31 // mesh-explorer -m ../../data/escher-p3.mesh
32 // mesh-explorer -m ../../data/mobius-strip.mesh
33 
34 #include "mfem.hpp"
35 #include <fstream>
36 #include <limits>
37 #include <cstdlib>
38 
39 using namespace mfem;
40 using namespace std;
41 
42 // This tranformation can be applied to a mesh with the 't' menu option.
43 void transformation(const Vector &p, Vector &v)
44 {
45  // simple shear transformation
46  double s = 0.1;
47 
48  if (p.Size() == 3)
49  {
50  v(0) = p(0) + s*p(1) + s*p(2);
51  v(1) = p(1) + s*p(2) + s*p(0);
52  v(2) = p(2);
53  }
54  else if (p.Size() == 2)
55  {
56  v(0) = p(0) + s*p(1);
57  v(1) = p(1) + s*p(0);
58  }
59  else
60  {
61  v = p;
62  }
63 }
64 
65 // This function is used with the 'r' menu option, sub-option 'l' to refine a
66 // mesh locally in a region, defined by return values <= region_eps.
67 double region_eps = 1e-8;
68 double region(const Vector &p)
69 {
70  const double x = p(0), y = p(1);
71  // here we describe the region: (x <= 1/4) && (y >= 0) && (y <= 1)
72  return std::max(std::max(x - 0.25, -y), y - 1.0);
73 }
74 
75 Mesh *read_par_mesh(int np, const char *mesh_prefix)
76 {
77  Mesh *mesh;
78  Array<Mesh *> mesh_array;
79 
80  mesh_array.SetSize(np);
81  for (int p = 0; p < np; p++)
82  {
83  ostringstream fname;
84  fname << mesh_prefix << '.' << setfill('0') << setw(6) << p;
85  ifgzstream meshin(fname.str().c_str());
86  if (!meshin)
87  {
88  cerr << "Can not open mesh file: " << fname.str().c_str()
89  << '!' << endl;
90  for (p--; p >= 0; p--)
91  {
92  delete mesh_array[p];
93  }
94  return NULL;
95  }
96  mesh_array[p] = new Mesh(meshin, 1, 0);
97  // set element and boundary attributes to be the processor number + 1
98  if (1)
99  {
100  for (int i = 0; i < mesh_array[p]->GetNE(); i++)
101  {
102  mesh_array[p]->GetElement(i)->SetAttribute(p+1);
103  }
104  for (int i = 0; i < mesh_array[p]->GetNBE(); i++)
105  {
106  mesh_array[p]->GetBdrElement(i)->SetAttribute(p+1);
107  }
108  }
109  }
110  mesh = new Mesh(mesh_array, np);
111 
112  for (int p = 0; p < np; p++)
113  {
114  delete mesh_array[np-1-p];
115  }
116  mesh_array.DeleteAll();
117 
118  return mesh;
119 }
120 
121 // Given a 3D mesh, produce a 2D mesh consisting of its boundary elements.
123 {
124  // Determine mapping from vertex to boundary vertex
125  Array<int> v2v(mesh->GetNV());
126  v2v = -1;
127  for (int i = 0; i < mesh->GetNBE(); i++)
128  {
129  Element *el = mesh->GetBdrElement(i);
130  int *v = el->GetVertices();
131  int nv = el->GetNVertices();
132  for (int j = 0; j < nv; j++)
133  {
134  v2v[v[j]] = 0;
135  }
136  }
137  int nbvt = 0;
138  for (int i = 0; i < v2v.Size(); i++)
139  {
140  if (v2v[i] == 0)
141  {
142  v2v[i] = nbvt++;
143  }
144  }
145 
146  // Create a new mesh for the boundary
147  Mesh * bmesh = new Mesh(mesh->Dimension() - 1, nbvt, mesh->GetNBE(),
148  0, mesh->SpaceDimension());
149 
150  // Copy vertices to the boundary mesh
151  nbvt = 0;
152  for (int i = 0; i < v2v.Size(); i++)
153  {
154  if (v2v[i] >= 0)
155  {
156  double *c = mesh->GetVertex(i);
157  bmesh->AddVertex(c);
158  nbvt++;
159  }
160  }
161 
162  // Copy elements to the boundary mesh
163  int bv[4];
164  for (int i = 0; i < mesh->GetNBE(); i++)
165  {
166  Element *el = mesh->GetBdrElement(i);
167  int *v = el->GetVertices();
168  int nv = el->GetNVertices();
169 
170  for (int j = 0; j < nv; j++)
171  {
172  bv[j] = v2v[v[j]];
173  }
174 
175  switch (el->GetGeometryType())
176  {
177  case Geometry::SEGMENT:
178  bmesh->AddSegment(bv, el->GetAttribute());
179  break;
180  case Geometry::TRIANGLE:
181  bmesh->AddTriangle(bv, el->GetAttribute());
182  break;
183  case Geometry::SQUARE:
184  bmesh->AddQuad(bv, el->GetAttribute());
185  break;
186  default:
187  break; /// This should not happen
188  }
189 
190  }
191  bmesh->FinalizeTopology();
192 
193  // Copy GridFunction describing nodes if present
194  if (mesh->GetNodes())
195  {
196  FiniteElementSpace *fes = mesh->GetNodes()->FESpace();
197  const FiniteElementCollection *fec = fes->FEColl();
198  if (dynamic_cast<const H1_FECollection*>(fec))
199  {
200  FiniteElementCollection *fec_copy =
202  FiniteElementSpace *fes_copy =
203  new FiniteElementSpace(*fes, bmesh, fec_copy);
204  GridFunction *bdr_nodes = new GridFunction(fes_copy);
205  bdr_nodes->MakeOwner(fec_copy);
206 
207  bmesh->NewNodes(*bdr_nodes, true);
208 
209  Array<int> vdofs;
210  Array<int> bvdofs;
211  Vector v;
212  for (int i=0; i<mesh->GetNBE(); i++)
213  {
214  fes->GetBdrElementVDofs(i, vdofs);
215  mesh->GetNodes()->GetSubVector(vdofs, v);
216 
217  fes_copy->GetElementVDofs(i, bvdofs);
218  bdr_nodes->SetSubVector(bvdofs, v);
219  }
220  }
221  else
222  {
223  cout << "\nDiscontinuous nodes not yet supported" << endl;
224  }
225  }
226 
227  return bmesh;
228 }
229 
230 int main (int argc, char *argv[])
231 {
232  int np = 0;
233  const char *mesh_file = "../../data/beam-hex.mesh";
234  bool refine = true;
235 
236  OptionsParser args(argc, argv);
237  args.AddOption(&mesh_file, "-m", "--mesh",
238  "Mesh file to visualize.");
239  args.AddOption(&np, "-np", "--num-proc",
240  "Load mesh from multiple processors.");
241  args.AddOption(&refine, "-ref", "--refinement", "-no-ref", "--no-refinement",
242  "Prepare the mesh for refinement or not.");
243  args.Parse();
244  if (!args.Good())
245  {
246  if (!args.Help())
247  {
248  args.PrintError(cout);
249  cout << endl;
250  }
251  cout << "Visualize and manipulate a serial mesh:\n"
252  << " mesh-explorer -m <mesh_file>\n"
253  << "Visualize and manipulate a parallel mesh:\n"
254  << " mesh-explorer -np <#proc> -m <mesh_prefix>\n" << endl
255  << "All Options:\n";
256  args.PrintHelp(cout);
257  return 1;
258  }
259  args.PrintOptions(cout);
260 
261  Mesh *mesh;
262  Mesh *bdr_mesh = NULL;
263  if (np <= 0)
264  {
265  mesh = new Mesh(mesh_file, 1, refine);
266  }
267  else
268  {
269  mesh = read_par_mesh(np, mesh_file);
270  if (mesh == NULL)
271  {
272  return 3;
273  }
274  }
275  int dim = mesh->Dimension();
276  int sdim = mesh->SpaceDimension();
277 
278  FiniteElementCollection *bdr_attr_fec = NULL;
279  FiniteElementCollection *attr_fec;
280  if (dim == 2)
281  {
282  attr_fec = new Const2DFECollection;
283  }
284  else
285  {
286  bdr_attr_fec = new Const2DFECollection;
287  attr_fec = new Const3DFECollection;
288  }
289 
290  int print_char = 1;
291  while (1)
292  {
293  if (print_char)
294  {
295  cout << endl;
296  mesh->PrintCharacteristics();
297  cout << "boundary attribs :";
298  for (int i = 0; i < mesh->bdr_attributes.Size(); i++)
299  {
300  cout << ' ' << mesh->bdr_attributes[i];
301  }
302  cout << '\n' << "material attribs :";
303  for (int i = 0; i < mesh->attributes.Size(); i++)
304  {
305  cout << ' ' << mesh->attributes[i];
306  }
307  cout << endl;
308  cout << "mesh curvature : ";
309  if (mesh->GetNodalFESpace() != NULL)
310  {
311  cout << mesh->GetNodalFESpace()->FEColl()->Name() << endl;
312  }
313  else
314  {
315  cout << "NONE" << endl;
316  }
317  }
318  print_char = 0;
319  cout << endl;
320  cout << "What would you like to do?\n"
321  "r) Refine\n"
322  "c) Change curvature\n"
323  "s) Scale\n"
324  "t) Transform\n"
325  "j) Jitter\n"
326  "v) View\n"
327  "m) View materials\n"
328  "b) View boundary\n"
329  "e) View elements\n"
330  "h) View element sizes, h\n"
331  "k) View element ratios, kappa\n"
332  "x) Print sub-element stats\n"
333  "f) Find physical point in reference space\n"
334  "p) Generate a partitioning\n"
335  "o) Reorder elements\n"
336  "S) Save in MFEM format\n"
337  "V) Save in VTK format (only linear and quadratic meshes)\n"
338  "q) Quit\n"
339 #ifdef MFEM_USE_ZLIB
340  "Z) Save in MFEM format with compression\n"
341 #endif
342  "--> " << flush;
343  char mk;
344  cin >> mk;
345  if (!cin) { break; }
346 
347  if (mk == 'q')
348  {
349  break;
350  }
351 
352  if (mk == 'r')
353  {
354  cout <<
355  "Choose type of refinement:\n"
356  "s) standard refinement with Mesh::UniformRefinement()\n"
357  "b) Mesh::UniformRefinement() (bisection for tet meshes)\n"
358  "u) uniform refinement with a factor\n"
359  "g) non-uniform refinement (Gauss-Lobatto) with a factor\n"
360  "l) refine locally using the region() function\n"
361  "--> " << flush;
362  char sk;
363  cin >> sk;
364  switch (sk)
365  {
366  case 's':
367  mesh->UniformRefinement();
368  // Make sure tet-only meshes are marked for local refinement.
369  mesh->Finalize(true);
370  break;
371  case 'b':
372  mesh->UniformRefinement(1); // ref_algo = 1
373  break;
374  case 'u':
375  case 'g':
376  {
377  cout << "enter refinement factor --> " << flush;
378  int ref_factor;
379  cin >> ref_factor;
380  if (ref_factor <= 1 || ref_factor > 32) { break; }
381  int ref_type = (sk == 'u') ? BasisType::ClosedUniform :
383  Mesh *rmesh = new Mesh(mesh, ref_factor, ref_type);
384  delete mesh;
385  mesh = rmesh;
386  break;
387  }
388  case 'l':
389  {
390  Vector pt;
391  Array<int> marked_elements;
392  for (int i = 0; i < mesh->GetNE(); i++)
393  {
394  // check all nodes of the element
396  mesh->GetElementTransformation(i, &T);
397  for (int j = 0; j < T.GetPointMat().Width(); j++)
398  {
399  T.GetPointMat().GetColumnReference(j, pt);
400  if (region(pt) <= region_eps)
401  {
402  marked_elements.Append(i);
403  break;
404  }
405  }
406  }
407  mesh->GeneralRefinement(marked_elements);
408  break;
409  }
410  }
411  print_char = 1;
412  }
413 
414  if (mk == 'c')
415  {
416  int p;
417  cout << "enter new order for mesh curvature --> " << flush;
418  cin >> p;
419  mesh->SetCurvature(p > 0 ? p : -p, p <= 0);
420  print_char = 1;
421  }
422 
423  if (mk == 's')
424  {
425  double factor;
426  cout << "scaling factor ---> " << flush;
427  cin >> factor;
428 
429  GridFunction *nodes = mesh->GetNodes();
430  if (nodes == NULL)
431  {
432  for (int i = 0; i < mesh->GetNV(); i++)
433  {
434  double *v = mesh->GetVertex(i);
435  v[0] *= factor;
436  v[1] *= factor;
437  if (dim == 3)
438  {
439  v[2] *= factor;
440  }
441  }
442  }
443  else
444  {
445  *nodes *= factor;
446  }
447 
448  print_char = 1;
449  }
450 
451  if (mk == 't')
452  {
453  mesh->Transform(transformation);
454  print_char = 1;
455  }
456 
457  if (mk == 'j')
458  {
459  double jitter;
460  cout << "jitter factor ---> " << flush;
461  cin >> jitter;
462 
463  GridFunction *nodes = mesh->GetNodes();
464 
465  if (nodes == NULL)
466  {
467  cerr << "The mesh should have nodes, introduce curvature first!\n";
468  }
469  else
470  {
471  FiniteElementSpace *fespace = nodes->FESpace();
472 
473  GridFunction rdm(fespace);
474  rdm.Randomize();
475  rdm -= 0.5; // shift to random values in [-0.5,0.5]
476  rdm *= jitter;
477 
478  // compute minimal local mesh size
479  Vector h0(fespace->GetNDofs());
480  h0 = infinity();
481  {
482  Array<int> dofs;
483  for (int i = 0; i < fespace->GetNE(); i++)
484  {
485  fespace->GetElementDofs(i, dofs);
486  for (int j = 0; j < dofs.Size(); j++)
487  {
488  h0(dofs[j]) = std::min(h0(dofs[j]), mesh->GetElementSize(i));
489  }
490  }
491  }
492 
493  // scale the random values to be of order of the local mesh size
494  for (int i = 0; i < fespace->GetNDofs(); i++)
495  {
496  for (int d = 0; d < dim; d++)
497  {
498  rdm(fespace->DofToVDof(i,d)) *= h0(i);
499  }
500  }
501 
502  char move_bdr = 'n';
503  cout << "move boundary nodes? [y/n] ---> " << flush;
504  cin >> move_bdr;
505 
506  // don't perturb the boundary
507  if (move_bdr == 'n')
508  {
509  Array<int> vdofs;
510  for (int i = 0; i < fespace->GetNBE(); i++)
511  {
512  fespace->GetBdrElementVDofs(i, vdofs);
513  for (int j = 0; j < vdofs.Size(); j++)
514  {
515  rdm(vdofs[j]) = 0.0;
516  }
517  }
518  }
519 
520  *nodes += rdm;
521  }
522 
523  print_char = 1;
524  }
525 
526  if (mk == 'x')
527  {
528  int sd, nz = 0;
529  DenseMatrix J(dim);
530  double min_det_J, max_det_J, min_det_J_z, max_det_J_z;
531  double min_kappa, max_kappa, max_ratio_det_J_z;
532  min_det_J = min_kappa = infinity();
533  max_det_J = max_kappa = max_ratio_det_J_z = -infinity();
534  cout << "subdivision factor ---> " << flush;
535  cin >> sd;
536  Array<int> bad_elems_by_geom(Geometry::NumGeom);
537  bad_elems_by_geom = 0;
538  for (int i = 0; i < mesh->GetNE(); i++)
539  {
542 
543  RefinedGeometry *RefG = GlobGeometryRefiner.Refine(geom, sd, 1);
544  IntegrationRule &ir = RefG->RefPts;
545 
546  min_det_J_z = infinity();
547  max_det_J_z = -infinity();
548  for (int j = 0; j < ir.GetNPoints(); j++)
549  {
550  T->SetIntPoint(&ir.IntPoint(j));
551  Geometries.JacToPerfJac(geom, T->Jacobian(), J);
552 
553  double det_J = J.Det();
554  double kappa =
555  J.CalcSingularvalue(0) / J.CalcSingularvalue(dim-1);
556 
557  min_det_J_z = fmin(min_det_J_z, det_J);
558  max_det_J_z = fmax(max_det_J_z, det_J);
559 
560  min_kappa = fmin(min_kappa, kappa);
561  max_kappa = fmax(max_kappa, kappa);
562  }
563  max_ratio_det_J_z =
564  fmax(max_ratio_det_J_z, max_det_J_z/min_det_J_z);
565  min_det_J = fmin(min_det_J, min_det_J_z);
566  max_det_J = fmax(max_det_J, max_det_J_z);
567  if (min_det_J_z <= 0.0)
568  {
569  nz++;
570  bad_elems_by_geom[geom]++;
571  }
572  }
573  cout << "\nbad elements = " << nz;
574  if (nz)
575  {
576  cout << " -- ";
577  Mesh::PrintElementsByGeometry(dim, bad_elems_by_geom, cout);
578  }
579  cout << "\nmin det(J) = " << min_det_J
580  << "\nmax det(J) = " << max_det_J
581  << "\nglobal ratio = " << max_det_J/min_det_J
582  << "\nmax el ratio = " << max_ratio_det_J_z
583  << "\nmin kappa = " << min_kappa
584  << "\nmax kappa = " << max_kappa << endl;
585  }
586 
587  if (mk == 'f')
588  {
589  DenseMatrix point_mat(sdim,1);
590  cout << "\npoint in physical space ---> " << flush;
591  for (int i = 0; i < sdim; i++)
592  {
593  cin >> point_mat(i,0);
594  }
595  Array<int> elem_ids;
597 
598  // physical -> reference space
599  mesh->FindPoints(point_mat, elem_ids, ips);
600 
601  cout << "point in reference space:";
602  if (elem_ids[0] == -1)
603  {
604  cout << " NOT FOUND!\n";
605  }
606  else
607  {
608  cout << " element " << elem_ids[0] << ", ip =";
609  cout << " " << ips[0].x;
610  if (sdim > 1)
611  {
612  cout << " " << ips[0].y;
613  if (sdim > 2)
614  {
615  cout << " " << ips[0].z;
616  }
617  }
618  cout << endl;
619  }
620  }
621 
622  if (mk == 'o')
623  {
624  cout << "What type of reordering?\n"
625  "h) Hilbert spatial sort\n"
626  //"g) Gecko edge-product minimization\n" // TODO future
627  "--> " << flush;
628  char rk;
629  cin >> rk;
630 
631  Array<int> ordering;
632  if (rk == 'h')
633  {
634  mesh->GetHilbertElementOrdering(ordering);
635  mesh->ReorderElements(ordering);
636  }
637  }
638 
639  // These are the cases that open a new GLVis window
640  if (mk == 'm' || mk == 'b' || mk == 'e' || mk == 'v' || mk == 'h' ||
641  mk == 'k' || mk == 'p')
642  {
643  Array<int> bdr_part;
644  Array<int> part(mesh->GetNE());
645  FiniteElementSpace *bdr_attr_fespace = NULL;
646  FiniteElementSpace *attr_fespace =
647  new FiniteElementSpace(mesh, attr_fec);
648  GridFunction bdr_attr;
649  GridFunction attr(attr_fespace);
650 
651  if (mk == 'm')
652  {
653  for (int i = 0; i < mesh->GetNE(); i++)
654  {
655  part[i] = (attr(i) = mesh->GetAttribute(i)) - 1;
656  }
657  }
658 
659  if (mk == 'b')
660  {
661  if (dim == 3)
662  {
663  delete bdr_mesh;
664  bdr_mesh = skin_mesh(mesh);
665  bdr_attr_fespace =
666  new FiniteElementSpace(bdr_mesh, bdr_attr_fec);
667  bdr_part.SetSize(bdr_mesh->GetNE());
668  bdr_attr.SetSpace(bdr_attr_fespace);
669  for (int i = 0; i < bdr_mesh->GetNE(); i++)
670  {
671  bdr_part[i] = (bdr_attr(i) = bdr_mesh->GetAttribute(i)) - 1;
672  }
673  }
674  else
675  {
676  attr = 1.0;
677  }
678  }
679 
680  if (mk == 'v')
681  {
682  attr = 1.0;
683  }
684 
685  if (mk == 'e')
686  {
687  Array<int> coloring;
688  srand(time(0));
689  double a = double(rand()) / (double(RAND_MAX) + 1.);
690  int el0 = (int)floor(a * mesh->GetNE());
691  cout << "Generating coloring starting with element " << el0+1
692  << " / " << mesh->GetNE() << endl;
693  mesh->GetElementColoring(coloring, el0);
694  for (int i = 0; i < coloring.Size(); i++)
695  {
696  attr(i) = coloring[i];
697  }
698  cout << "Number of colors: " << attr.Max() + 1 << endl;
699  for (int i = 0; i < mesh->GetNE(); i++)
700  {
701  // part[i] = i; // checkerboard element coloring
702  attr(i) = part[i] = i; // coloring by element number
703  }
704  }
705 
706  if (mk == 'h')
707  {
708  DenseMatrix J(dim);
709  double h_min, h_max;
710  h_min = infinity();
711  h_max = -h_min;
712  for (int i = 0; i < mesh->GetNE(); i++)
713  {
714  int geom = mesh->GetElementBaseGeometry(i);
716  T->SetIntPoint(&Geometries.GetCenter(geom));
717  Geometries.JacToPerfJac(geom, T->Jacobian(), J);
718 
719  attr(i) = J.Det();
720  if (attr(i) < 0.0)
721  {
722  attr(i) = -pow(-attr(i), 1.0/double(dim));
723  }
724  else
725  {
726  attr(i) = pow(attr(i), 1.0/double(dim));
727  }
728  h_min = min(h_min, attr(i));
729  h_max = max(h_max, attr(i));
730  }
731  cout << "h_min = " << h_min << ", h_max = " << h_max << endl;
732  }
733 
734  if (mk == 'k')
735  {
736  DenseMatrix J(dim);
737  for (int i = 0; i < mesh->GetNE(); i++)
738  {
739  int geom = mesh->GetElementBaseGeometry(i);
741  T->SetIntPoint(&Geometries.GetCenter(geom));
742  Geometries.JacToPerfJac(geom, T->Jacobian(), J);
743  attr(i) = J.CalcSingularvalue(0) / J.CalcSingularvalue(dim-1);
744  }
745  }
746 
747  if (mk == 'p')
748  {
749  int *partitioning = NULL, np;
750  cout << "What type of partitioning?\n"
751  "c) Cartesian\n"
752  "s) Simple 1D split of the element sequence\n"
753  "0) METIS_PartGraphRecursive (sorted neighbor lists)\n"
754  "1) METIS_PartGraphKway (sorted neighbor lists)"
755  " (default)\n"
756  "2) METIS_PartGraphVKway (sorted neighbor lists)\n"
757  "3) METIS_PartGraphRecursive\n"
758  "4) METIS_PartGraphKway\n"
759  "5) METIS_PartGraphVKway\n"
760  "--> " << flush;
761  char pk;
762  cin >> pk;
763  if (pk == 'c')
764  {
765  int nxyz[3];
766  cout << "Enter nx: " << flush;
767  cin >> nxyz[0]; np = nxyz[0];
768  if (mesh->Dimension() > 1)
769  {
770  cout << "Enter ny: " << flush;
771  cin >> nxyz[1]; np *= nxyz[1];
772  if (mesh->Dimension() > 2)
773  {
774  cout << "Enter nz: " << flush;
775  cin >> nxyz[2]; np *= nxyz[2];
776  }
777  }
778  partitioning = mesh->CartesianPartitioning(nxyz);
779  }
780  else if (pk == 's')
781  {
782  cout << "Enter number of processors: " << flush;
783  cin >> np;
784 
785  partitioning = new int[mesh->GetNE()];
786  for (int i = 0; i < mesh->GetNE(); i++)
787  {
788  partitioning[i] = i * np / mesh->GetNE();
789  }
790  }
791  else
792  {
793  int part_method = pk - '0';
794  if (part_method < 0 || part_method > 5)
795  {
796  continue;
797  }
798  cout << "Enter number of processors: " << flush;
799  cin >> np;
800  partitioning = mesh->GeneratePartitioning(np, part_method);
801  }
802  if (partitioning)
803  {
804  const char part_file[] = "partitioning.txt";
805  ofstream opart(part_file);
806  opart << "number_of_elements " << mesh->GetNE() << '\n'
807  << "number_of_processors " << np << '\n';
808  for (int i = 0; i < mesh->GetNE(); i++)
809  {
810  opart << partitioning[i] << '\n';
811  }
812  cout << "Partitioning file: " << part_file << endl;
813 
814  Array<int> proc_el(np);
815  proc_el = 0;
816  for (int i = 0; i < mesh->GetNE(); i++)
817  {
818  proc_el[partitioning[i]]++;
819  }
820  int min_el = proc_el[0], max_el = proc_el[0];
821  for (int i = 1; i < np; i++)
822  {
823  if (min_el > proc_el[i])
824  {
825  min_el = proc_el[i];
826  }
827  if (max_el < proc_el[i])
828  {
829  max_el = proc_el[i];
830  }
831  }
832  cout << "Partitioning stats:\n"
833  << " "
834  << setw(12) << "minimum"
835  << setw(12) << "average"
836  << setw(12) << "maximum"
837  << setw(12) << "total" << '\n';
838  cout << " elements "
839  << setw(12) << min_el
840  << setw(12) << double(mesh->GetNE())/np
841  << setw(12) << max_el
842  << setw(12) << mesh->GetNE() << endl;
843  }
844  else
845  {
846  continue;
847  }
848 
849  for (int i = 0; i < mesh->GetNE(); i++)
850  {
851  attr(i) = part[i] = partitioning[i];
852  }
853  delete [] partitioning;
854  }
855 
856  char vishost[] = "localhost";
857  int visport = 19916;
858  socketstream sol_sock(vishost, visport);
859  if (sol_sock.is_open())
860  {
861  sol_sock.precision(14);
862  if (sdim == 2)
863  {
864  sol_sock << "fem2d_gf_data_keys\n";
865  if (mk != 'p')
866  {
867  mesh->Print(sol_sock);
868  }
869  else
870  {
871  // NURBS meshes do not support PrintWithPartitioning
872  if (mesh->NURBSext)
873  {
874  mesh->Print(sol_sock);
875  for (int i = 0; i < mesh->GetNE(); i++)
876  {
877  attr(i) = part[i];
878  }
879  }
880  else
881  {
882  mesh->PrintWithPartitioning(part, sol_sock, 1);
883  }
884  }
885  attr.Save(sol_sock);
886  sol_sock << "RjlmAb***********";
887  if (mk == 'v')
888  {
889  sol_sock << "e";
890  }
891  else
892  {
893  sol_sock << "\n";
894  }
895  }
896  else
897  {
898  sol_sock << "fem3d_gf_data_keys\n";
899  if (mk == 'v' || mk == 'h' || mk == 'k')
900  {
901  mesh->Print(sol_sock);
902  }
903  else if (mk == 'b')
904  {
905  bdr_mesh->Print(sol_sock);
906  bdr_attr.Save(sol_sock);
907  sol_sock << "mcaaA";
908  // Switch to a discrete color scale
909  sol_sock << "pppppp" << "pppppp" << "pppppp";
910  }
911  else
912  {
913  // NURBS meshes do not support PrintWithPartitioning
914  if (mesh->NURBSext)
915  {
916  mesh->Print(sol_sock);
917  for (int i = 0; i < mesh->GetNE(); i++)
918  {
919  attr(i) = part[i];
920  }
921  }
922  else
923  {
924  mesh->PrintWithPartitioning(part, sol_sock);
925  }
926  }
927  if (mk != 'b')
928  {
929  attr.Save(sol_sock);
930  sol_sock << "maaA";
931  if (mk == 'v')
932  {
933  sol_sock << "aa";
934  }
935  else
936  {
937  sol_sock << "\n";
938  }
939  }
940  }
941  sol_sock << flush;
942  }
943  else
944  {
945  cout << "Unable to connect to "
946  << vishost << ':' << visport << endl;
947  }
948  delete attr_fespace;
949  delete bdr_attr_fespace;
950  }
951 
952  if (mk == 'S')
953  {
954  const char mesh_file[] = "mesh-explorer.mesh";
955  ofstream omesh(mesh_file);
956  omesh.precision(14);
957  mesh->Print(omesh);
958  cout << "New mesh file: " << mesh_file << endl;
959  }
960 
961  if (mk == 'V')
962  {
963  const char mesh_file[] = "mesh-explorer.vtk";
964  ofstream omesh(mesh_file);
965  omesh.precision(14);
966  mesh->PrintVTK(omesh);
967  cout << "New VTK mesh file: " << mesh_file << endl;
968  }
969 
970 #ifdef MFEM_USE_ZLIB
971  if (mk == 'Z')
972  {
973  const char mesh_file[] = "mesh-explorer.mesh.gz";
974  ofgzstream omesh(mesh_file, "zwb9");
975  omesh.precision(14);
976  mesh->Print(omesh);
977  cout << "New mesh file: " << mesh_file << endl;
978  }
979 #endif
980 
981  }
982 
983  delete bdr_attr_fec;
984  delete attr_fec;
985  delete bdr_mesh;
986  delete mesh;
987  return 0;
988 }
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:245
bool Help() const
Definition: optparser.hpp:126
Geometry::Type GetGeometryType() const
Definition: element.hpp:52
int Size() const
Logical size of the array.
Definition: array.hpp:124
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
Definition: vector.cpp:498
virtual void Print(std::ostream &out=mfem::out) const
Definition: mesh.hpp:1188
int * CartesianPartitioning(int nxyz[])
Definition: mesh.cpp:5312
int GetNDofs() const
Returns number of degrees of freedom.
Definition: fespace.hpp:378
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
DenseMatrix & GetPointMat()
Read and write access to the underlying point matrix describing the transformation.
Definition: eltrans.hpp:321
const double * GetVertex(int i) const
Return pointer to vertex i&#39;s coordinates.
Definition: mesh.hpp:750
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:27
int DofToVDof(int dof, int vd, int ndofs=-1) const
Definition: fespace.cpp:143
Mesh * read_par_mesh(int np, const char *mesh_prefix)
virtual void GetVertices(Array< int > &v) const =0
Returns element&#39;s vertices.
static void PrintElementsByGeometry(int dim, const Array< int > &num_elems_by_geom, std::ostream &out)
Auxiliary method used by PrintCharacteristics().
Definition: mesh.cpp:219
static const int NumGeom
Definition: geom.hpp:41
void JacToPerfJac(int GeomType, const DenseMatrix &J, DenseMatrix &PJ) const
Definition: geom.cpp:737
int * GeneratePartitioning(int nparts, int part_method=1)
Definition: mesh.cpp:5351
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:696
void PrintWithPartitioning(int *partitioning, std::ostream &out, int elem_attr=0) const
Prints the mesh with boundary elements given by the boundary of the subdomains, so that the boundary ...
Definition: mesh.cpp:9166
void MakeOwner(FiniteElementCollection *_fec)
Make the GridFunction the owner of fec and fes.
Definition: gridfunc.hpp:112
const Geometry::Type geom
Definition: ex1.cpp:40
double Det() const
Definition: densemat.cpp:449
void NewNodes(GridFunction &nodes, bool make_owner=false)
Replace the internal node GridFunction with the given GridFunction.
Definition: mesh.cpp:6227
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:63
void SetIntPoint(const IntegrationPoint *ip)
Definition: eltrans.hpp:56
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
int Size() const
Returns the size of the vector.
Definition: vector.hpp:157
void Transform(void(*f)(const Vector &, Vector &))
Definition: mesh.cpp:9864
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:693
void PrintHelp(std::ostream &out) const
Definition: optparser.cpp:378
void Randomize(int seed=0)
Set random values in the vector.
Definition: vector.cpp:692
const IntegrationPoint & GetCenter(int GeomType)
Return the center of the given Geometry::Type, GeomType.
Definition: geom.hpp:70
int main(int argc, char *argv[])
Definition: ex1.cpp:62
Geometry::Type GetElementBaseGeometry(int i) const
Definition: mesh.hpp:790
void DeleteAll()
Delete whole array.
Definition: array.hpp:802
Mesh * skin_mesh(Mesh *mesh)
void AddVertex(const double *)
Definition: mesh.cpp:1155
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:2731
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:408
Geometry Geometries
Definition: fe.cpp:12638
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:248
int GetNBE() const
Returns number of boundary elements in the mesh.
Definition: fespace.hpp:417
void AddSegment(const int *vi, int attr=1)
Definition: mesh.cpp:1166
int Append(const T &el)
Append element to array, resize if necessary.
Definition: array.hpp:707
void transformation(const Vector &p, Vector &v)
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:7982
double region_eps
const DenseMatrix & Jacobian()
Return the Jacobian matrix of the transformation at the currently set IntegrationPoint, using the method SetIntPoint().
Definition: eltrans.hpp:71
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:3924
virtual void GetElementDofs(int i, Array< int > &dofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: fespace.cpp:1540
GeometryRefiner GlobGeometryRefiner
Definition: geom.cpp:1395
int GetAttribute() const
Return element&#39;s attribute.
Definition: element.hpp:55
void GetHilbertElementOrdering(Array< int > &ordering)
Definition: mesh.cpp:1552
IntegrationRule RefPts
Definition: geom.hpp:239
RefinedGeometry * Refine(Geometry::Type Geom, int Times, int ETimes=1)
Definition: geom.cpp:925
Nodes: x_i = i/(n-1), i=0,...,n-1.
Definition: fe.hpp:37
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:442
int Dimension() const
Definition: mesh.hpp:744
double GetElementSize(int i, int type=0)
Get the size of the i-th element relative to the perfect reference element.
Definition: mesh.cpp:76
int SpaceDimension() const
Definition: mesh.hpp:745
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:189
void PrintVTK(std::ostream &out)
Definition: mesh.cpp:8446
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:87
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)
Definition: optparser.hpp:76
void FinalizeTopology(bool generate_bdr=true)
Finalize the construction of the secondary topology (connectivity) data of a Mesh.
Definition: mesh.cpp:2242
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:635
void AddQuad(const int *vi, int attr=1)
Definition: mesh.cpp:1181
void GetColumnReference(int c, Vector &col)
Definition: densemat.hpp:270
static FiniteElementCollection * New(const char *name)
Factory method: return a newly allocated FiniteElementCollection according to the given name...
Definition: fe_coll.cpp:47
double a
Definition: lissajous.cpp:41
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:191
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
Definition: mesh.hpp:690
void PrintError(std::ostream &out) const
Definition: optparser.cpp:335
virtual const char * Name() const
Definition: fe_coll.hpp:53
void AddTriangle(const int *vi, int attr=1)
Definition: mesh.cpp:1176
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
Definition: mesh.cpp:338
double Max() const
Returns the maximal element of the vector.
Definition: vector.cpp:799
const FiniteElementSpace * GetNodalFESpace() const
Definition: mesh.cpp:3919
virtual void Finalize(bool refine=false, bool fix_orientation=false)
Finalize the construction of a general Mesh.
Definition: mesh.cpp:2316
void ReorderElements(const Array< int > &ordering, bool reorder_vertices=true)
Definition: mesh.cpp:1604
int dim
Definition: ex24.cpp:43
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
double CalcSingularvalue(const int i) const
Return the i-th singular value (decreasing order) of NxN matrix, N=1,2,3.
Definition: densemat.cpp:1224
double infinity()
Define a shortcut for std::numeric_limits&lt;double&gt;::infinity()
Definition: vector.hpp:42
virtual int GetNVertices() const =0
virtual void SetSpace(FiniteElementSpace *f)
Associate a new FiniteElementSpace with the GridFunction.
Definition: gridfunc.cpp:180
virtual int FindPoints(DenseMatrix &point_mat, Array< int > &elem_ids, Array< IntegrationPoint > &ips, bool warn=true, InverseElementTransformation *inv_trans=NULL)
Find the ids of the elements that contain the given points, and their corresponding reference coordin...
Definition: mesh.cpp:10133
double kappa
Definition: ex3.cpp:54
Vector data type.
Definition: vector.hpp:48
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:395
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:6203
void GetElementColoring(Array< int > &colors, int el0=0)
Definition: mesh.cpp:9091
void GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i&#39;th boundary element.
Definition: fespace.cpp:178
void GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
Definition: mesh.cpp:7577
double region(const Vector &p)
Abstract data type element.
Definition: element.hpp:28
int GetAttribute(int i) const
Return the attribute of element i.
Definition: mesh.hpp:1001
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Definition: mesh.hpp:187
const Element * GetBdrElement(int i) const
Definition: mesh.hpp:779
void PrintCharacteristics(Vector *Vh=NULL, Vector *Vk=NULL, std::ostream &out=mfem::out)
Compute and print mesh characteristics such as number of vertices, number of elements, number of boundary elements, minimal and maximal element sizes, minimal and maximal element aspect ratios, etc.
Definition: mesh.cpp:233
bool Good() const
Definition: optparser.hpp:125