MFEM  v4.5.2
Finite element discretization library
mondrian.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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 // Mondrian Miniapp: Convert an image to an AMR mesh
14 // -------------------------------------------------
15 //
16 // This miniapp is a specialized version of the Shaper miniapp that converts an
17 // input image to an AMR mesh. It allows the fast approximate meshing of any
18 // domain for which there is an image.
19 //
20 // The input to image should be in 8-bit grayscale PGM format. You can use a
21 // number of image manipulation tools, such as GIMP (gimp.org) and ImageMagick's
22 // convert utility (imagemagick.org/script/convert.php) to convert your image to
23 // this format as a pre-processing step, e.g.:
24 //
25 // /usr/bin/convert australia.svg -compress none -depth 8 australia.pgm
26 //
27 // Compile with: make mondrian
28 //
29 // Sample runs: mondrian -i australia.pgm
30 // mondrian -i australia.pgm -m ../../data/inline-tri.mesh
31 // mondrian -i australia.pgm -m ../../data/disc-nurbs.mesh
32 // mondrian -i australia.pgm -sd 3 -a -ncl -1
33 
34 #include "mfem.hpp"
35 #include <fstream>
36 #include <iostream>
37 
38 using namespace mfem;
39 using namespace std;
40 
41 // Simple class to parse portable graymap format (PGM) image files, see
42 // http://netpbm.sourceforge.net/doc/pgm.html
43 class ParsePGM
44 {
45 public:
46  ParsePGM(const char *filename);
47  ~ParsePGM();
48 
49  int Height() const { return N; }
50  int Width() const { return M; }
51 
52  int operator()(int i, int j) const
53  { return int((pgm8) ? pgm8[M*i+j] : pgm16[M*i+j]); }
54 
55 private:
56  int M, N;
57  int depth;
58 
59  char *pgm8;
60  unsigned short int *pgm16;
61 
62  void ReadMagicNumber(istream &in);
63  void ReadComments(istream &in);
64  void ReadDimensions(istream &in);
65  void ReadDepth(istream &in);
66  void ReadPGM(istream &in);
67 };
68 
69 // Given a point x, return its "material" specification defined by the grayscale
70 // pixel values from the pgm image using NC different colors.
71 int material(const ParsePGM &pgm, int NC,
72  Vector &x, Vector &xmin, Vector &xmax);
73 
74 int main(int argc, char *argv[])
75 {
76  const char *mesh_file = "../../data/inline-quad.mesh";
77  const char *img_file = "australia.pgm";
78  int sd = 2;
79  int nclimit = 1;
80  int ncolors = 3;
81  bool aniso = false;
82  bool visualization = 1;
83 
84  // Parse command line
85  OptionsParser args(argc, argv);
86  args.AddOption(&mesh_file, "-m", "--mesh",
87  "Input mesh file to shape materials in.");
88  args.AddOption(&img_file, "-i", "--img",
89  "Input image.");
90  args.AddOption(&sd, "-sd", "--sub-divisions",
91  "Number of element subdivisions for interface detection.");
92  args.AddOption(&nclimit, "-ncl", "--nc-limit",
93  "Level of hanging nodes allowed (-1 = unlimited).");
94  args.AddOption(&ncolors, "-nc", "--num-colors",
95  "Number of colors considered (1-256, based on binning).");
96  args.AddOption(&aniso, "-a", "--aniso", "-i", "--iso",
97  "Enable anisotropic refinement of quads and hexes.");
98  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
99  "--no-visualization",
100  "Enable or disable GLVis visualization.");
101  args.Parse();
102  if (!args.Good()) { args.PrintUsage(cout); return 1; }
103  args.PrintOptions(cout);
104 
105  // Read the image
106  ParsePGM pgm(img_file);
107 
108  // Read initial mesh, get dimensions and bounding box
109  Mesh mesh(mesh_file, 1, 1);
110  int dim = mesh.Dimension();
111  int sdim = mesh.SpaceDimension();
112  Vector xmin, xmax;
113  mesh.GetBoundingBox(xmin, xmax);
114 
115  // NURBS meshes don't support non-conforming refinement for now
116  if (mesh.NURBSext) { mesh.SetCurvature(2); }
117 
118  // Anisotropic refinement not supported for simplex meshes.
119  if (mesh.MeshGenerator() & 1) { aniso = false; }
120 
121  // Mesh attributes will be visualized as piece-wise constants
122  L2_FECollection attr_fec(0, dim);
123  FiniteElementSpace attr_fespace(&mesh, &attr_fec);
124  GridFunction attr(&attr_fespace);
125 
126  // GLVis server to visualize to
127  socketstream sol_sock;
128  if (visualization)
129  {
130  char vishost[] = "localhost";
131  int visport = 19916;
132  sol_sock.open(vishost, visport);
133  sol_sock.precision(8);
134  }
135 
136  // Shaping loop
137  for (int iter = 0; 1; iter++)
138  {
139  Array<Refinement> refs;
140  for (int i = 0; i < mesh.GetNE(); i++)
141  {
142  bool refine = false;
143 
144  // Sample materials in each element using "sd" sub-divisions
145  Vector pt;
146  Geometry::Type geom = mesh.GetElementBaseGeometry(i);
148  RefinedGeometry *RefG = GlobGeometryRefiner.Refine(geom, sd, 1);
149  IntegrationRule &ir = RefG->RefPts;
150 
151  // Refine any element where different materials are detected. A more
152  // sophisticated logic can be implemented here -- e.g. don't refine
153  // the interfaces between certain materials.
154  Array<int> mat(ir.GetNPoints());
155  double matsum = 0.0;
156  for (int j = 0; j < ir.GetNPoints(); j++)
157  {
158  T->Transform(ir.IntPoint(j), pt);
159  int m = material(pgm, 256/ncolors, pt, xmin, xmax);
160  mat[j] = m;
161  matsum += m;
162  if ((int)matsum != m*(j+1))
163  {
164  refine = true;
165  }
166  }
167 
168  // Set the element attribute as the "average". Other choices are
169  // possible here too, e.g. attr(i) = mat;
170  attr(i) = round(matsum/ir.GetNPoints());
171 
172  // Mark the element for refinement
173  if (refine)
174  {
175  int type = 7;
176  if (aniso)
177  {
178  // Determine the XYZ bitmask for anisotropic refinement.
179  int dx = 0, dy = 0, dz = 0;
180  const int s = sd+1;
181  if (dim == 2)
182  {
183  for (int jj = 0; jj <= sd; jj++)
184  for (int ii = 0; ii < sd; ii++)
185  {
186  dx += abs(mat[jj*s + ii+1] - mat[jj*s + ii]);
187  dy += abs(mat[(ii+1)*s + jj] - mat[ii*s + jj]);
188  }
189  }
190  else if (dim == 3)
191  {
192  for (int kk = 0; kk <= sd; kk++)
193  for (int jj = 0; jj <= sd; jj++)
194  for (int ii = 0; ii < sd; ii++)
195  {
196  dx += abs(mat[(kk*s + jj)*s + ii+1] - mat[(kk*s + jj)*s + ii]);
197  dy += abs(mat[(kk*s + ii+1)*s + jj] - mat[(kk*s + ii)*s + jj]);
198  dz += abs(mat[((ii+1)*s + jj)*s + kk] - mat[(ii*s + jj)*s + kk]);
199  }
200  }
201  type = 0;
202  const int tol = mat.Size() / 10;
203  if (dx > tol) { type |= 1; }
204  if (dy > tol) { type |= 2; }
205  if (dz > tol) { type |= 4; }
206  if (!type) { type = 7; } // because of tol
207  }
208 
209  refs.Append(Refinement(i, type));
210  }
211  }
212 
213  // Visualization
214  if (visualization)
215  {
216  sol_sock << "solution\n" << mesh << attr;
217  if (iter == 0 && sdim == 2)
218  {
219  sol_sock << "keys 'RjlmpppppppppppppA*************'\n";
220  }
221  if (iter == 0 && sdim == 3)
222  {
223  sol_sock << "keys 'YYYYYYYYYXXXXXXXmA********8888888pppttt";
224  if (dim == 3) { sol_sock << "iiM"; }
225  sol_sock << "'\n";
226  }
227  sol_sock << flush;
228  }
229 
230  // Ask the user if we should continue refining
231  cout << "Iteration " << iter+1 << ": mesh has " << mesh.GetNE() <<
232  " elements. \n";
233  if ((iter+1) % 3 == 0)
234  {
235  if (!visualization) { break; }
236  char yn;
237  cout << "Continue shaping? --> ";
238  cin >> yn;
239  if (yn == 'n' || yn == 'q') { break; }
240  }
241 
242  // Perform refinement, update spaces and grid functions
243  mesh.GeneralRefinement(refs, -1, nclimit);
244  attr_fespace.Update();
245  attr.Update();
246  }
247 
248  // Set element attributes in the mesh object before saving
249  for (int i = 0; i < mesh.GetNE(); i++)
250  {
251  mesh.SetAttribute(i, static_cast<int>(attr(i)));
252  }
253  mesh.SetAttributes();
254 
255  // Save the final mesh
256  ofstream mesh_ofs("mondrian.mesh");
257  mesh_ofs.precision(8);
258  mesh.Print(mesh_ofs);
259 }
260 
261 ParsePGM::ParsePGM(const char *filename)
262  : M(-1), N(-1), depth(-1), pgm8(NULL), pgm16(NULL)
263 {
264  ifstream in(filename);
265  if (!in)
266  {
267  // Abort with an error message
268  MFEM_ABORT("Image file not found: " << filename << '\n');
269  }
270 
271  ReadMagicNumber(in);
272  ReadDimensions(in);
273  ReadDepth(in);
274  ReadPGM(in);
275 
276  in.close();
277 }
278 
279 ParsePGM::~ParsePGM()
280 {
281  if (pgm8 != NULL) { delete [] pgm8; }
282  if (pgm16 != NULL) { delete [] pgm16; }
283 }
284 
285 void ParsePGM::ReadMagicNumber(istream &in)
286 {
287  char c;
288  int p;
289  in >> c >> p; // Read magic number which should be P2 or P5
290  MFEM_VERIFY(c == 'P' && (p == 2 || p == 5),
291  "Invalid PGM file! Unrecognized magic number\""
292  << c << p << "\".");
293  ReadComments(in);
294 }
295 
296 void ParsePGM::ReadComments(istream &in)
297 {
298  string buf;
299  in >> std::ws; // absorb any white space
300  while (in.peek() == '#')
301  {
302  std::getline(in,buf);
303  }
304  in >> std::ws; // absorb any white space
305 }
306 
307 void ParsePGM::ReadDimensions(istream &in)
308 {
309  in >> M;
310  ReadComments(in);
311  in >> N;
312  ReadComments(in);
313 }
314 
315 void ParsePGM::ReadDepth(istream &in)
316 {
317  in >> depth;
318  ReadComments(in);
319 }
320 
321 void ParsePGM::ReadPGM(istream &in)
322 {
323  if (depth < 16)
324  {
325  pgm8 = new char[M*N];
326  }
327  else
328  {
329  pgm16 = new unsigned short int[M*N];
330  }
331 
332  if (pgm8)
333  {
334  for (int i=0; i<M*N; i++)
335  {
336  in >> pgm8[i];
337  }
338  }
339  else
340  {
341  for (int i=0; i<M*N; i++)
342  {
343  in >> pgm16[i];
344  }
345  }
346 }
347 
348 int material(const ParsePGM &pgm, int NC, Vector &x, Vector &xmin, Vector &xmax)
349 {
350  // Rescaling to [0,1]^sdim
351  for (int i = 0; i < x.Size(); i++)
352  {
353  x(i) = (x(i)-xmin(i))/(xmax(i)-xmin(i));
354  }
355 
356  int M = pgm.Width();
357  int N = pgm.Height();
358 
359  int i = (int)(x(1)*N), j = (int)(x(0)*M);
360  if (i == N) { i = N-1; }
361  if (j == M) { j = M-1; }
362  i = N-1-i;
363 
364  return pgm(i,j)/NC+1;
365 }
const char vishost[]
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
virtual void Update(bool want_transform=true)
Reflect changes in the mesh: update number of DOFs, etc. Also, calculate GridFunction transformation ...
Definition: fespace.cpp:3336
int material(const ParsePGM &pgm, int NC, Vector &x, Vector &xmin, Vector &xmax)
Definition: mondrian.cpp:348
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
Geometry::Type GetElementBaseGeometry(int i) const
Definition: mesh.hpp:1108
int Dimension() const
Definition: mesh.hpp:1047
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
void GetBoundingBox(Vector &min, Vector &max, int ref=2)
Returns the minimum and maximum corners of the mesh bounding box.
Definition: mesh.cpp:129
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150
STL namespace.
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:250
GeometryRefiner GlobGeometryRefiner
Definition: geom.cpp:1773
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
virtual void SetAttributes()
Definition: mesh.cpp:1564
int Append(const T &el)
Append element &#39;el&#39; to array, resize if necessary.
Definition: array.hpp:756
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Set the curvature of the mesh nodes using the given polynomial degree.
Definition: mesh.cpp:5356
int MeshGenerator()
Get the mesh generator/type.
Definition: mesh.hpp:929
IntegrationRule RefPts
Definition: geom.hpp:314
RefinedGeometry * Refine(Geometry::Type Geom, int Times, int ETimes=1)
Definition: geom.cpp:1099
const int visport
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
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
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Definition: gridfunc.cpp:165
int SpaceDimension() const
Definition: mesh.hpp:1048
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:936
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:277
int main(int argc, char *argv[])
Definition: mondrian.cpp:74
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
Definition: mesh.cpp:348
int dim
Definition: ex24.cpp:53
int open(const char hostname[], int port)
Open the socket stream on &#39;port&#39; at &#39;hostname&#39;.
void SetAttribute(int i, int attr)
Set the attribute of element i.
Definition: mesh.hpp:1547
Vector data type.
Definition: vector.hpp:60
virtual void Print(std::ostream &os=mfem::out) const
Definition: mesh.hpp:1749
RefCoord s[3]
void GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
Definition: mesh.cpp:9467
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:320