MFEM  v4.1.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
shaper.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 // Shaper Miniapp: Resolve material interfaces by mesh refinement
14 // --------------------------------------------------------------
15 //
16 // This miniapp performs multiple levels of adaptive mesh refinement to resolve
17 // the interfaces between different "materials" in the mesh, as specified by the
18 // given material() function. It can be used as a simple initial mesh generator,
19 // for example in the case when the interface is too complex to describe without
20 // local refinement. Both conforming and non-conforming refinements are supported.
21 //
22 // Two additional versions of this miniapp can be found in the miniapps/toys
23 // directory: Mandel uses the Shaper algorithm for fractal visualization, while
24 // Mondrian convert an image to an AMR mesh suitable for MFEM computations.
25 //
26 // Compile with: make shaper
27 //
28 // Sample runs: shaper
29 // shaper -m ../../data/inline-tri.mesh
30 // shaper -m ../../data/inline-hex.mesh
31 // shaper -m ../../data/inline-tet.mesh
32 // shaper -m ../../data/amr-quad.mesh
33 // shaper -m ../../data/beam-quad.mesh -a -ncl -1 -sd 4
34 // shaper -m ../../data/ball-nurbs.mesh
35 // shaper -m ../../data/mobius-strip.mesh
36 // shaper -m ../../data/square-disc-surf.mesh
37 // shaper -m ../../data/star-q3.mesh -sd 2 -ncl -1
38 // shaper -m ../../data/fichera-amr.mesh -a -ncl -1
39 
40 #include "mfem.hpp"
41 #include <fstream>
42 #include <iostream>
43 
44 using namespace mfem;
45 using namespace std;
46 
47 // Given a point x, return its material id as an integer. The ids should be
48 // positive. If the point is exactly on the interface, return 0.
49 //
50 // This particular implementation, rescales the mesh to [-1,1]^sdim given the
51 // xmin/xmax bounding box, and shapes in a simple annulus/shell with respect to
52 // the rescaled coordinates.
53 int material(Vector &x, Vector &xmin, Vector &xmax)
54 {
55  static double p = 2.0;
56 
57  // Rescaling to [-1,1]^sdim
58  for (int i = 0; i < x.Size(); i++)
59  {
60  x(i) = (2*x(i)-xmin(i)-xmax(i))/(xmax(i)-xmin(i));
61  }
62 
63  // A simple annulus/shell
64  if (x.Normlp(p) > 0.4 && x.Normlp(p) < 0.6) { return 1; }
65  if (x.Normlp(p) < 0.4 || x.Normlp(p) > 0.6) { return 2; }
66  return 0;
67 }
68 
69 int main(int argc, char *argv[])
70 {
71  int sd = 2;
72  int nclimit = 1;
73  const char *mesh_file = "../../data/inline-quad.mesh";
74  bool aniso = false;
75 
76  // Parse command line
77  OptionsParser args(argc, argv);
78  args.AddOption(&mesh_file, "-m", "--mesh",
79  "Input mesh file to shape materials in.");
80  args.AddOption(&sd, "-sd", "--sub-divisions",
81  "Number of element subdivisions for interface detection.");
82  args.AddOption(&nclimit, "-ncl", "--nc-limit",
83  "Level of hanging nodes allowed (-1 = unlimited).");
84  args.AddOption(&aniso, "-a", "--aniso", "-i", "--iso",
85  "Enable anisotropic refinement of quads and hexes.");
86  args.Parse();
87  if (!args.Good()) { args.PrintUsage(cout); return 1; }
88  args.PrintOptions(cout);
89 
90  // Read initial mesh, get dimensions and bounding box
91  Mesh mesh(mesh_file, 1, 1);
92  int dim = mesh.Dimension();
93  int sdim = mesh.SpaceDimension();
94  Vector xmin, xmax;
95  mesh.GetBoundingBox(xmin, xmax);
96 
97  // NURBS meshes don't support non-conforming refinement for now
98  if (mesh.NURBSext) { mesh.SetCurvature(2); }
99 
100  // Anisotropic refinement not supported for simplex meshes.
101  if (mesh.MeshGenerator() & 1) { aniso = false; }
102 
103  // Mesh attributes will be visualized as piece-wise constants
104  L2_FECollection attr_fec(0, dim);
105  FiniteElementSpace attr_fespace(&mesh, &attr_fec);
106  GridFunction attr(&attr_fespace);
107 
108  // GLVis server to visualize to
109  char vishost[] = "localhost";
110  int visport = 19916;
111  socketstream sol_sock(vishost, visport);
112  sol_sock.precision(8);
113 
114  // Shaping loop
115  for (int iter = 0; 1; iter++)
116  {
117  Array<Refinement> refs;
118  for (int i = 0; i < mesh.GetNE(); i++)
119  {
120  bool refine = false;
121 
122  // Sample materials in each element using "sd" sub-divisions
123  Vector pt;
126  RefinedGeometry *RefG = GlobGeometryRefiner.Refine(geom, sd, 1);
127  IntegrationRule &ir = RefG->RefPts;
128 
129  // Refine any element where different materials are detected. A more
130  // sophisticated logic can be implemented here -- e.g. don't refine
131  // the interfaces between certain materials.
132  Array<int> mat(ir.GetNPoints());
133  double matsum = 0.0;
134  for (int j = 0; j < ir.GetNPoints(); j++)
135  {
136  T->Transform(ir.IntPoint(j), pt);
137  int m = material(pt, xmin, xmax);
138  mat[j] = m;
139  matsum += m;
140  if ((int)matsum != m*(j+1))
141  {
142  refine = true;
143  }
144  }
145 
146  // Set the element attribute as the "average". Other choices are
147  // possible here too, e.g. attr(i) = mat;
148  attr(i) = round(matsum/ir.GetNPoints());
149 
150  // Mark the element for refinement
151  if (refine)
152  {
153  int type = 7;
154  if (aniso)
155  {
156  // Determine the XYZ bitmask for anisotropic refinement.
157  int dx = 0, dy = 0, dz = 0;
158  const int s = sd+1;
159  if (dim == 2)
160  {
161  for (int j = 0; j <= sd; j++)
162  for (int i = 0; i < sd; i++)
163  {
164  dx += abs(mat[j*s + i+1] - mat[j*s + i]);
165  dy += abs(mat[(i+1)*s + j] - mat[i*s + j]);
166  }
167  }
168  else if (dim == 3)
169  {
170  for (int k = 0; k <= sd; k++)
171  for (int j = 0; j <= sd; j++)
172  for (int i = 0; i < sd; i++)
173  {
174  dx += abs(mat[(k*s + j)*s + i+1] - mat[(k*s + j)*s + i]);
175  dy += abs(mat[(k*s + i+1)*s + j] - mat[(k*s + i)*s + j]);
176  dz += abs(mat[((i+1)*s + j)*s + k] - mat[(i*s + j)*s + k]);
177  }
178  }
179  type = 0;
180  const int tol = mat.Size() / 10;
181  if (dx > tol) { type |= 1; }
182  if (dy > tol) { type |= 2; }
183  if (dz > tol) { type |= 4; }
184  if (!type) { type = 7; } // because of tol
185  }
186 
187  refs.Append(Refinement(i, type));
188  }
189  }
190 
191  // Visualization
192  sol_sock << "solution\n" << mesh << attr;
193  if (iter == 0 && sdim == 2)
194  {
195  sol_sock << "keys 'RjlmpppppppppppppA*************'\n";
196  }
197  if (iter == 0 && sdim == 3)
198  {
199  sol_sock << "keys 'YYYYYYYYYXXXXXXXmA********8888888pppttt";
200  if (dim == 3) { sol_sock << "iiM"; }
201  sol_sock << "'\n";
202  }
203  sol_sock << flush;
204 
205  // Ask the user if we should continue refining
206  char yn;
207  cout << "Mesh has " << mesh.GetNE() << " elements. \n"
208  << "Continue shaping? --> ";
209  cin >> yn;
210  if (yn == 'n' || yn == 'q') { break; }
211 
212  // Perform refinement, update spaces and grid functions
213  mesh.GeneralRefinement(refs, -1, nclimit);
214  attr_fespace.Update();
215  attr.Update();
216  }
217 
218  // Set element attributes in the mesh object before saving
219  for (int i = 0; i < mesh.GetNE(); i++)
220  {
221  mesh.SetAttribute(i, attr(i));
222  }
223  mesh.SetAttributes();
224 
225  // Save the final mesh
226  ofstream mesh_ofs("shaper.mesh");
227  mesh_ofs.precision(8);
228  mesh.Print(mesh_ofs);
229 }
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
double Normlp(double p) const
Returns the l_p norm of the vector.
Definition: vector.cpp:749
virtual void SetAttributes()
Definition: mesh.cpp:1107
int Append(const T &el)
Append element to array, resize if necessary.
Definition: array.hpp:707
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
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
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