MFEM v4.8.0
Finite element discretization library
Loading...
Searching...
No Matches
mandel.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2025, 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
31using namespace mfem;
32using 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.
39int 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 real_t col = p(0), row = p(1);
48 {
49 int width = 1080, height = 1080;
50 col *= width;
51 row *= height;
52 real_t c_re = (col - width/2)*4.0/width;
53 real_t c_im = (row - height/2)*4.0/width;
54 real_t x = 0, y = 0;
55 int iteration = 0, maxit = 10000;
56 while (x*x+y*y <= 4 && iteration < maxit)
57 {
58 real_t 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
74int 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 int visport = 19916;
82
83 // Parse command line
84 OptionsParser args(argc, argv);
85 args.AddOption(&mesh_file, "-m", "--mesh",
86 "Input mesh file to shape materials in.");
87 args.AddOption(&sd, "-sd", "--sub-divisions",
88 "Number of element subdivisions for interface detection.");
89 args.AddOption(&nclimit, "-ncl", "--nc-limit",
90 "Level of hanging nodes allowed (-1 = unlimited).");
91 args.AddOption(&aniso, "-a", "--aniso", "-i", "--iso",
92 "Enable anisotropic refinement of quads and hexes.");
93 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
94 "--no-visualization",
95 "Enable or disable GLVis visualization.");
96 args.AddOption(&visport, "-p", "--send-port", "Socket for GLVis.");
97 args.Parse();
98 if (!args.Good()) { args.PrintUsage(cout); return 1; }
99 args.PrintOptions(cout);
100
101 // Read initial mesh, get dimensions and bounding box
102 Mesh mesh(mesh_file, 1, 1);
103 int dim = mesh.Dimension();
104 int sdim = mesh.SpaceDimension();
105 Vector xmin, xmax;
106 mesh.GetBoundingBox(xmin, xmax);
107
108 // Increase the mesh resolution
109 for (int l = 0; l < 3; l++) { mesh.UniformRefinement(); }
110
111 // NURBS meshes don't support non-conforming refinement for now
112 if (mesh.NURBSext) { mesh.SetCurvature(2); }
113
114 // Anisotropic refinement not supported for simplex meshes.
115 if (mesh.MeshGenerator() & 1) { aniso = false; }
116
117 // Mesh attributes will be visualized as piece-wise constants
118 L2_FECollection attr_fec(0, dim);
119 FiniteElementSpace attr_fespace(&mesh, &attr_fec);
120 GridFunction attr(&attr_fespace);
121
122 // GLVis server to visualize to
123 socketstream sol_sock;
124 if (visualization)
125 {
126 char vishost[] = "localhost";
127 sol_sock.open(vishost, visport);
128 sol_sock.precision(8);
129 }
130
131 // Shaping loop
132 for (int iter = 0; 1; iter++)
133 {
135 for (int e = 0; e < mesh.GetNE(); e++)
136 {
137 bool refine = false;
138
139 // Sample materials in each element using "sd" sub-divisions
140 Vector pt;
143 RefinedGeometry *RefG = GlobGeometryRefiner.Refine(geom, sd, 1);
144 IntegrationRule &ir = RefG->RefPts;
145
146 // Refine any element where different materials are detected. A more
147 // sophisticated logic can be implemented here -- e.g. don't refine
148 // the interfaces between certain materials.
149 Array<int> mat(ir.GetNPoints());
150 real_t matsum = 0.0;
151 for (int j = 0; j < ir.GetNPoints(); j++)
152 {
153 T->Transform(ir.IntPoint(j), pt);
154 int m = material(pt, xmin, xmax);
155 mat[j] = m;
156 matsum += m;
157 if ((int)matsum != m*(j+1))
158 {
159 refine = true;
160 }
161 }
162
163 // Set the element attribute as the "average". Other choices are
164 // possible here too, e.g. attr(e) = mat;
165 attr(e) = round(matsum/ir.GetNPoints());
166
167 // Mark the element for refinement
168 if (refine)
169 {
170 int type = 7;
171 if (aniso)
172 {
173 // Determine the XYZ bitmask for anisotropic refinement.
174 int dx = 0, dy = 0, dz = 0;
175 const int s = sd+1;
176 if (dim == 2)
177 {
178 for (int j = 0; j <= sd; j++)
179 for (int i = 0; i < sd; i++)
180 {
181 dx += abs(mat[j*s + i+1] - mat[j*s + i]);
182 dy += abs(mat[(i+1)*s + j] - mat[i*s + j]);
183 }
184 }
185 else if (dim == 3)
186 {
187 for (int k = 0; k <= sd; k++)
188 for (int j = 0; j <= sd; j++)
189 for (int i = 0; i < sd; i++)
190 {
191 dx += abs(mat[(k*s + j)*s + i+1] - mat[(k*s + j)*s + i]);
192 dy += abs(mat[(k*s + i+1)*s + j] - mat[(k*s + i)*s + j]);
193 dz += abs(mat[((i+1)*s + j)*s + k] - mat[(i*s + j)*s + k]);
194 }
195 }
196 type = 0;
197 const int tol = mat.Size() / 10;
198 if (dx > tol) { type |= 1; }
199 if (dy > tol) { type |= 2; }
200 if (dz > tol) { type |= 4; }
201 if (!type) { type = 7; } // because of tol
202 }
203
204 refs.Append(Refinement(e, type));
205 }
206 }
207
208 // Visualization
209 if (visualization)
210 {
211 sol_sock << "solution\n" << mesh << attr;
212 if (iter == 0 && sdim == 2)
213 {
214 sol_sock << "keys 'RjlppppppppppppppA*************'\n";
215 }
216 if (iter == 0 && sdim == 3)
217 {
218 sol_sock << "keys 'YYYYYYYYYXXXXXXXmA********8888888pppttt";
219 if (dim == 3) { sol_sock << "iiM"; }
220 sol_sock << "'\n";
221 }
222 sol_sock << flush;
223 }
224
225 // Ask the user if we should continue refining
226 cout << "Iteration " << iter+1 << ": mesh has " << mesh.GetNE() <<
227 " elements. \n";
228 if ((iter+1) % 4 == 0)
229 {
230 if (!visualization) { break; }
231 char yn;
232 cout << "Continue shaping? --> ";
233 cin >> yn;
234 if (yn == 'n' || yn == 'q') { break; }
235 }
236
237 // Perform refinement, update spaces and grid functions
238 mesh.GeneralRefinement(refs, -1, nclimit);
239 attr_fespace.Update();
240 attr.Update();
241 }
242
243 // Set element attributes in the mesh object before saving
244 for (int i = 0; i < mesh.GetNE(); i++)
245 {
246 mesh.SetAttribute(i, static_cast<int>(attr(i)));
247 }
248 mesh.SetAttributes();
249
250 // Save the final mesh
251 ofstream mesh_ofs("mandel.mesh");
252 mesh_ofs.precision(8);
253 mesh.Print(mesh_ofs);
254}
int Size() const
Return the logical size of the array.
Definition array.hpp:147
int Append(const T &el)
Append element 'el' to array, resize if necessary.
Definition array.hpp:830
virtual void Transform(const IntegrationPoint &, Vector &)=0
Transform integration point from reference coordinates to physical coordinates and store them in the ...
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:244
virtual void Update(bool want_transform=true)
Reflect changes in the mesh: update number of DOFs, etc. Also, calculate GridFunction transformation ...
Definition fespace.cpp:4157
RefinedGeometry * Refine(Geometry::Type Geom, int Times, int ETimes=1)
Definition geom.cpp:1136
Class for grid function - Vector with associated FE space.
Definition gridfunc.hpp:31
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Definition gridfunc.cpp:167
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:100
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition intrules.hpp:256
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition intrules.hpp:259
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition fe_coll.hpp:346
Mesh data type.
Definition mesh.hpp:64
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition mesh.hpp:298
void GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
Definition mesh.cpp:10883
void SetAttribute(int i, int attr)
Set the attribute of element i.
Definition mesh.cpp:7629
virtual void Print(std::ostream &os=mfem::out, const std::string &comments="") const
Definition mesh.hpp:2433
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1282
void GetBoundingBox(Vector &min, Vector &max, int ref=2)
Returns the minimum and maximum corners of the mesh bounding box.
Definition mesh.cpp:138
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1216
void GetElementTransformation(int i, IsoparametricTransformation *ElTr) const
Builds the transformation defining the i-th element in ElTr. ElTr must be allocated in advance and wi...
Definition mesh.cpp:357
int SpaceDimension() const
Dimension of the physical space containing the mesh.
Definition mesh.hpp:1219
int MeshGenerator() const
Get the mesh generator/type.
Definition mesh.hpp:1243
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:6484
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:11295
Geometry::Type GetElementBaseGeometry(int i) const
Definition mesh.hpp:1455
virtual void SetAttributes()
Determine the sets of unique attribute values in domain and boundary elements.
Definition mesh.cpp:1819
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void PrintUsage(std::ostream &out) const
Print the usage message.
void PrintOptions(std::ostream &out) const
Print the options.
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 'var' to receive the value. Enable/disable tags are used to set the bool...
Definition optparser.hpp:82
bool Good() const
Return true if the command line options were parsed successfully.
IntegrationRule RefPts
Definition geom.hpp:321
Vector data type.
Definition vector.hpp:82
int open(const char hostname[], int port)
Open the socket stream on 'port' at 'hostname'.
int dim
Definition ex24.cpp:53
int main()
int material(Vector &p, Vector &pmin, Vector &pmax)
Definition mandel.cpp:39
GeometryRefiner GlobGeometryRefiner
Definition geom.cpp:2014
float real_t
Definition config.hpp:43
const char vishost[]
real_t p(const Vector &x, real_t t)