MFEM  v4.1.0 Finite element discretization library
mandel.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
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 &x, Vector &xmin, Vector &xmax)
40 {
41  // Rescaling to [0,1]^sdim
42  for (int i = 0; i < x.Size(); i++)
43  {
44  x(i) = (x(i)-xmin(i))/(xmax(i)-xmin(i));
45  }
46  x(0) -= 0.1;
47  double col = x(0), row = x(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);
85  "Input mesh file to shape materials in.");
87  "Number of element subdivisions for interface detection.");
89  "Level of hanging nodes allowed (-1 = unlimited).");
90  args.AddOption(&aniso, "-a", "--aniso", "-i", "--iso",
91  "Enable anisotropic refinement of quads and hexes.");
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 i = 0; i < mesh.GetNE(); i++)
135  {
136  bool refine = false;
137
138  // Sample materials in each element using "sd" sub-divisions
139  Vector pt;
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(i) = mat;
164  attr(i) = 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(i, 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, 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 }
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:245
virtual void Print(std::ostream &out=mfem::out) const
Definition: mesh.hpp:1188
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:27
virtual void Update(bool want_transform=true)
Definition: fespace.cpp:2057
const Geometry::Type geom
Definition: ex1.cpp:40
void GetBoundingBox(Vector &min, Vector &max, int ref=2)
Returns the minimum and maximum corners of the mesh bounding box.
Definition: mesh.cpp:120
int Size() const
Returns the size of the vector.
Definition: vector.hpp:157
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:693
int main(int argc, char *argv[])
Definition: ex1.cpp:62
Geometry::Type GetElementBaseGeometry(int i) const
Definition: mesh.hpp:790
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:248
int material(Vector &x, Vector &xmin, Vector &xmax)
Definition: shaper.cpp:53
virtual void SetAttributes()
Definition: mesh.cpp:1107
int Append(const T &el)
Append element to array, resize if necessary.
Definition: array.hpp:707
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:7982
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:3924
int MeshGenerator()
Get the mesh generator/type.
Definition: mesh.hpp:686
GeometryRefiner GlobGeometryRefiner
Definition: geom.cpp:1395
IntegrationRule RefPts
Definition: geom.hpp:239
RefinedGeometry * Refine(Geometry::Type Geom, int Times, int ETimes=1)
Definition: geom.cpp:925
int Dimension() const
Definition: mesh.hpp:744
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
int SpaceDimension() const
Definition: mesh.hpp:745
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
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Definition: gridfunc.cpp:152
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:191
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
Definition: mesh.cpp:338
int dim
Definition: ex24.cpp:43
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
int open(const char hostname[], int port)
void SetAttribute(int i, int attr)
Set the attribute of element i.
Definition: mesh.hpp:1004
Vector data type.
Definition: vector.hpp:48
virtual void Transform(const IntegrationPoint &, Vector &)=0
void GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
Definition: mesh.cpp:7577
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:143
bool Good() const
Definition: optparser.hpp:125