MFEM  v4.6.0
Finite element discretization library
mandel.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 // Mandel Miniapp: Fractal visualization with AMR
14 // ----------------------------------------------
15 //
16 // This miniapp is a specialized version of the Shaper miniapp to the Mandelbrot
17 // set. It provides a light-hearted example of AMR and GLVis integration.
18 //
19 // Compile with: make mandel
20 //
21 // Sample runs: mandel
22 // mandel -m ../../data/inline-tri.mesh
23 // mandel -m ../../data/star-mixed-p2.mesh -a -ncl -1 -sd 4
24 // mandel -m ../../data/klein-bottle.mesh
25 // mandel -m ../../data/inline-hex.mesh
26 
27 #include "mfem.hpp"
28 #include <fstream>
29 #include <iostream>
30 
31 using namespace mfem;
32 using namespace std;
33 
34 // Given a point x, return its material id as an integer. The ids should be
35 // positive. If the point is exactly on the interface, return 0.
36 //
37 // In this particular miniapp, the material value is based on the number of
38 // iterations for the point from the definition of the Mandelbrot set.
39 int material(Vector &p, Vector &pmin, Vector &pmax)
40 {
41  // Rescaling to [0,1]^sdim
42  for (int i = 0; i < p.Size(); i++)
43  {
44  p(i) = (p(i)-pmin(i))/(pmax(i)-pmin(i));
45  }
46  p(0) -= 0.1;
47  double col = p(0), row = p(1);
48  {
49  int width = 1080, height = 1080;
50  col *= width;
51  row *= height;
52  double c_re = (col - width/2)*4.0/width;
53  double c_im = (row - height/2)*4.0/width;
54  double x = 0, y = 0;
55  int iteration = 0, maxit = 10000;
56  while (x*x+y*y <= 4 && iteration < maxit)
57  {
58  double x_new = x*x - y*y + c_re;
59  y = 2*x*y + c_im;
60  x = x_new;
61  iteration++;
62  }
63  if (iteration < maxit)
64  {
65  return iteration%10+2;
66  }
67  else
68  {
69  return 1;
70  }
71  }
72 }
73 
74 int main(int argc, char *argv[])
75 {
76  const char *mesh_file = "../../data/inline-quad.mesh";
77  int sd = 2;
78  int nclimit = 1;
79  bool aniso = false;
80  bool visualization = 1;
81 
82  // Parse command line
83  OptionsParser args(argc, argv);
84  args.AddOption(&mesh_file, "-m", "--mesh",
85  "Input mesh file to shape materials in.");
86  args.AddOption(&sd, "-sd", "--sub-divisions",
87  "Number of element subdivisions for interface detection.");
88  args.AddOption(&nclimit, "-ncl", "--nc-limit",
89  "Level of hanging nodes allowed (-1 = unlimited).");
90  args.AddOption(&aniso, "-a", "--aniso", "-i", "--iso",
91  "Enable anisotropic refinement of quads and hexes.");
92  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
93  "--no-visualization",
94  "Enable or disable GLVis visualization.");
95  args.Parse();
96  if (!args.Good()) { args.PrintUsage(cout); return 1; }
97  args.PrintOptions(cout);
98 
99  // Read initial mesh, get dimensions and bounding box
100  Mesh mesh(mesh_file, 1, 1);
101  int dim = mesh.Dimension();
102  int sdim = mesh.SpaceDimension();
103  Vector xmin, xmax;
104  mesh.GetBoundingBox(xmin, xmax);
105 
106  // Increase the mesh resolution
107  for (int l = 0; l < 3; l++) { mesh.UniformRefinement(); }
108 
109  // NURBS meshes don't support non-conforming refinement for now
110  if (mesh.NURBSext) { mesh.SetCurvature(2); }
111 
112  // Anisotropic refinement not supported for simplex meshes.
113  if (mesh.MeshGenerator() & 1) { aniso = false; }
114 
115  // Mesh attributes will be visualized as piece-wise constants
116  L2_FECollection attr_fec(0, dim);
117  FiniteElementSpace attr_fespace(&mesh, &attr_fec);
118  GridFunction attr(&attr_fespace);
119 
120  // GLVis server to visualize to
121  socketstream sol_sock;
122  if (visualization)
123  {
124  char vishost[] = "localhost";
125  int visport = 19916;
126  sol_sock.open(vishost, visport);
127  sol_sock.precision(8);
128  }
129 
130  // Shaping loop
131  for (int iter = 0; 1; iter++)
132  {
133  Array<Refinement> refs;
134  for (int e = 0; e < mesh.GetNE(); e++)
135  {
136  bool refine = false;
137 
138  // Sample materials in each element using "sd" sub-divisions
139  Vector pt;
140  Geometry::Type geom = mesh.GetElementBaseGeometry(e);
142  RefinedGeometry *RefG = GlobGeometryRefiner.Refine(geom, sd, 1);
143  IntegrationRule &ir = RefG->RefPts;
144 
145  // Refine any element where different materials are detected. A more
146  // sophisticated logic can be implemented here -- e.g. don't refine
147  // the interfaces between certain materials.
148  Array<int> mat(ir.GetNPoints());
149  double matsum = 0.0;
150  for (int j = 0; j < ir.GetNPoints(); j++)
151  {
152  T->Transform(ir.IntPoint(j), pt);
153  int m = material(pt, xmin, xmax);
154  mat[j] = m;
155  matsum += m;
156  if ((int)matsum != m*(j+1))
157  {
158  refine = true;
159  }
160  }
161 
162  // Set the element attribute as the "average". Other choices are
163  // possible here too, e.g. attr(e) = mat;
164  attr(e) = round(matsum/ir.GetNPoints());
165 
166  // Mark the element for refinement
167  if (refine)
168  {
169  int type = 7;
170  if (aniso)
171  {
172  // Determine the XYZ bitmask for anisotropic refinement.
173  int dx = 0, dy = 0, dz = 0;
174  const int s = sd+1;
175  if (dim == 2)
176  {
177  for (int j = 0; j <= sd; j++)
178  for (int i = 0; i < sd; i++)
179  {
180  dx += abs(mat[j*s + i+1] - mat[j*s + i]);
181  dy += abs(mat[(i+1)*s + j] - mat[i*s + j]);
182  }
183  }
184  else if (dim == 3)
185  {
186  for (int k = 0; k <= sd; k++)
187  for (int j = 0; j <= sd; j++)
188  for (int i = 0; i < sd; i++)
189  {
190  dx += abs(mat[(k*s + j)*s + i+1] - mat[(k*s + j)*s + i]);
191  dy += abs(mat[(k*s + i+1)*s + j] - mat[(k*s + i)*s + j]);
192  dz += abs(mat[((i+1)*s + j)*s + k] - mat[(i*s + j)*s + k]);
193  }
194  }
195  type = 0;
196  const int tol = mat.Size() / 10;
197  if (dx > tol) { type |= 1; }
198  if (dy > tol) { type |= 2; }
199  if (dz > tol) { type |= 4; }
200  if (!type) { type = 7; } // because of tol
201  }
202 
203  refs.Append(Refinement(e, type));
204  }
205  }
206 
207  // Visualization
208  if (visualization)
209  {
210  sol_sock << "solution\n" << mesh << attr;
211  if (iter == 0 && sdim == 2)
212  {
213  sol_sock << "keys 'RjlppppppppppppppA*************'\n";
214  }
215  if (iter == 0 && sdim == 3)
216  {
217  sol_sock << "keys 'YYYYYYYYYXXXXXXXmA********8888888pppttt";
218  if (dim == 3) { sol_sock << "iiM"; }
219  sol_sock << "'\n";
220  }
221  sol_sock << flush;
222  }
223 
224  // Ask the user if we should continue refining
225  cout << "Iteration " << iter+1 << ": mesh has " << mesh.GetNE() <<
226  " elements. \n";
227  if ((iter+1) % 4 == 0)
228  {
229  if (!visualization) { break; }
230  char yn;
231  cout << "Continue shaping? --> ";
232  cin >> yn;
233  if (yn == 'n' || yn == 'q') { break; }
234  }
235 
236  // Perform refinement, update spaces and grid functions
237  mesh.GeneralRefinement(refs, -1, nclimit);
238  attr_fespace.Update();
239  attr.Update();
240  }
241 
242  // Set element attributes in the mesh object before saving
243  for (int i = 0; i < mesh.GetNE(); i++)
244  {
245  mesh.SetAttribute(i, static_cast<int>(attr(i)));
246  }
247  mesh.SetAttributes();
248 
249  // Save the final mesh
250  ofstream mesh_ofs("mandel.mesh");
251  mesh_ofs.precision(8);
252  mesh.Print(mesh_ofs);
253 }
const char vishost[]
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:253
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:96
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:3402
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:331
Geometry::Type GetElementBaseGeometry(int i) const
Definition: mesh.hpp:1242
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1020
int main(int argc, char *argv[])
Definition: mandel.cpp:74
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:462
int material(Vector &p, Vector &pmin, Vector &pmax)
Definition: mandel.cpp:39
void GetBoundingBox(Vector &min, Vector &max, int ref=2)
Returns the minimum and maximum corners of the mesh bounding box.
Definition: mesh.cpp:136
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:159
STL namespace.
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:256
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()
Determine the sets of unique attribute values in domain and boundary elements.
Definition: mesh.cpp:1572
int Append(const T &el)
Append element &#39;el&#39; to array, resize if necessary.
Definition: array.hpp:759
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:10232
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:5635
int MeshGenerator()
Get the mesh generator/type.
Definition: mesh.hpp:1047
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:219
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:164
int SpaceDimension() const
Dimension of the physical space containing the mesh.
Definition: mesh.hpp:1023
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:1086
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:277
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
Definition: mesh.cpp:355
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:1193
Vector data type.
Definition: vector.hpp:58
virtual void Print(std::ostream &os=mfem::out) const
Definition: mesh.hpp:2011
RefCoord s[3]
void GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
Definition: mesh.cpp:9820
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:327