MFEM  v3.3.2
Finite element discretization library
shaper.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
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 // Compile with: make shaper
23 //
24 // Sample runs: shaper
25 // shaper -m ../../data/inline-tri.mesh
26 // shaper -m ../../data/inline-hex.mesh
27 // shaper -m ../../data/inline-tet.mesh
28 // shaper -m ../../data/amr-quad.mesh
29 // shaper -m ../../data/beam-quad.mesh -a -ncl -1 -sd 4
30 // shaper -m ../../data/ball-nurbs.mesh
31 // shaper -m ../../data/mobius-strip.mesh
32 // shaper -m ../../data/square-disc-surf.mesh
33 // shaper -m ../../data/star-q3.mesh -sd 2 -ncl -1
34 // shaper -m ../../data/fichera-amr.mesh -a -ncl -1
35 
36 #include "mfem.hpp"
37 #include <fstream>
38 #include <iostream>
39 
40 using namespace mfem;
41 using namespace std;
42 
43 // Given a point x, return its material id as an integer. The ids should be
44 // positive. If the point is exactly on the interface, return 0.
45 //
46 // This particular implementation, rescales the mesh to [-1,1]^sdim given the
47 // xmin/xmax bounding box, and shapes in a simple annulus/shell with respect to
48 // the rescaled coordinates.
49 int material(Vector &x, Vector &xmin, Vector &xmax)
50 {
51  static double p = 2.0;
52 
53  // Rescaling to [-1,1]^sdim
54  for (int i = 0; i < x.Size(); i++)
55  {
56  x(i) = (2*x(i)-xmin(i)-xmax(i))/(xmax(i)-xmin(i));
57  }
58 
59  // A simple annulus/shell
60  if (x.Normlp(p) > 0.4 && x.Normlp(p) < 0.6) { return 1; }
61  if (x.Normlp(p) < 0.4 || x.Normlp(p) > 0.6) { return 2; }
62  return 0;
63 }
64 
65 int main(int argc, char *argv[])
66 {
67  int sd = 2;
68  int nclimit = 1;
69  const char *mesh_file = "../../data/inline-quad.mesh";
70  bool aniso = false;
71 
72  // Parse command line
73  OptionsParser args(argc, argv);
74  args.AddOption(&mesh_file, "-m", "--mesh",
75  "Input mesh file to shape materials in.");
76  args.AddOption(&sd, "-sd", "--sub-divisions",
77  "Number of element subdivisions for interface detection.");
78  args.AddOption(&nclimit, "-ncl", "--nc-limit",
79  "Level of hanging nodes allowed (-1 = unlimited).");
80  args.AddOption(&aniso, "-a", "--aniso", "-i", "--iso",
81  "Enable anisotropic refinement of quads and hexes.");
82  args.Parse();
83  if (!args.Good()) { args.PrintUsage(cout); return 1; }
84  args.PrintOptions(cout);
85 
86  // Read initial mesh, get dimensions and bounding box
87  Mesh mesh(mesh_file, 1, 1);
88  int dim = mesh.Dimension();
89  int sdim = mesh.SpaceDimension();
90  Vector xmin, xmax;
91  mesh.GetBoundingBox(xmin, xmax);
92 
93  // NURBS meshes don't support non-conforming refinement for now
94  if (mesh.NURBSext) { mesh.SetCurvature(2); }
95 
96  // Anisotropic refinement not supported for simplex meshes.
97  if (mesh.MeshGenerator() & 1) { aniso = false; }
98 
99  // Mesh attributes will be visualized as piece-wise constants
100  L2_FECollection attr_fec(0, dim);
101  FiniteElementSpace attr_fespace(&mesh, &attr_fec);
102  GridFunction attr(&attr_fespace);
103 
104  // GLVis server to visualize to
105  char vishost[] = "localhost";
106  int visport = 19916;
107  socketstream sol_sock(vishost, visport);
108  sol_sock.precision(8);
109 
110  // Shaping loop
111  for (int iter = 0; 1; iter++)
112  {
113  Array<Refinement> refs;
114  for (int i = 0; i < mesh.GetNE(); i++)
115  {
116  bool refine = false;
117 
118  // Sample materials in each element using "sd" sub-divisions
119  Vector pt;
120  int geom = mesh.GetElementBaseGeometry(i);
122  RefinedGeometry *RefG = GlobGeometryRefiner.Refine(geom, sd, 1);
123  IntegrationRule &ir = RefG->RefPts;
124 
125  // Refine any element where different materials are detected. A more
126  // sophisticated logic can be implemented here -- e.g. don't refine
127  // the interfaces between certain materials.
128  Array<int> mat(ir.GetNPoints());
129  double matsum = 0.0;
130  for (int j = 0; j < ir.GetNPoints(); j++)
131  {
132  T->Transform(ir.IntPoint(j), pt);
133  int m = material(pt, xmin, xmax);
134  mat[j] = m;
135  matsum += m;
136  if ((int)matsum != m*(j+1))
137  {
138  refine = true;
139  }
140  }
141 
142  // Set the element attribute as the "average". Other choices are
143  // possible here too, e.g. attr(i) = mat;
144  attr(i) = round(matsum/ir.GetNPoints());
145 
146  // Mark the element for refinement
147  if (refine)
148  {
149  int type = 7;
150  if (aniso)
151  {
152  // Determine the XYZ bitmask for anisotropic refinement.
153  int dx = 0, dy = 0, dz = 0;
154  const int s = sd+1;
155  if (dim == 2)
156  {
157  for (int j = 0; j <= sd; j++)
158  for (int i = 0; i < sd; i++)
159  {
160  dx += abs(mat[j*s + i+1] - mat[j*s + i]);
161  dy += abs(mat[(i+1)*s + j] - mat[i*s + j]);
162  }
163  }
164  else if (dim == 3)
165  {
166  for (int k = 0; k <= sd; k++)
167  for (int j = 0; j <= sd; j++)
168  for (int i = 0; i < sd; i++)
169  {
170  dx += abs(mat[(k*s + j)*s + i+1] - mat[(k*s + j)*s + i]);
171  dy += abs(mat[(k*s + i+1)*s + j] - mat[(k*s + i)*s + j]);
172  dz += abs(mat[((i+1)*s + j)*s + k] - mat[(i*s + j)*s + k]);
173  }
174  }
175  type = 0;
176  const int tol = mat.Size() / 10;
177  if (dx > tol) { type |= 1; }
178  if (dy > tol) { type |= 2; }
179  if (dz > tol) { type |= 4; }
180  if (!type) { type = 7; } // because of tol
181  }
182 
183  refs.Append(Refinement(i, type));
184  }
185  }
186 
187  // Visualization
188  sol_sock << "solution\n" << mesh << attr;
189  if (iter == 0 && sdim == 2)
190  {
191  sol_sock << "keys 'RjlmpppppppppppppA*************'\n";
192  }
193  if (iter == 0 && sdim == 3)
194  {
195  sol_sock << "keys 'YYYYYYYYYXXXXXXXmA********8888888pppttt";
196  if (dim == 3) { sol_sock << "iiM"; }
197  sol_sock << "'\n";
198  }
199  sol_sock << flush;
200 
201  // Ask the user if we should continue refining
202  char yn;
203  cout << "Mesh has " << mesh.GetNE() << " elements. \n"
204  << "Continue shaping? --> ";
205  cin >> yn;
206  if (yn == 'n' || yn == 'q') { break; }
207 
208  // Perform refinement, update spaces and grid functions
209  mesh.GeneralRefinement(refs, -1, nclimit);
210  attr_fespace.Update();
211  attr.Update();
212  }
213 
214  // Set element attributes in the mesh object before saving
215  for (int i = 0; i < mesh.GetNE(); i++)
216  {
217  mesh.SetAttribute(i, attr(i));
218  }
219  mesh.SetAttributes();
220 
221  // Save the final mesh
222  ofstream mesh_ofs("shaper.mesh");
223  mesh_ofs.precision(8);
224  mesh.Print(mesh_ofs);
225 }
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:222
virtual void Print(std::ostream &out=mfem::out) const
Definition: mesh.hpp:1019
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:83
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:27
virtual void Update(bool want_transform=true)
Definition: fespace.cpp:1494
int main(int argc, char *argv[])
Definition: shaper.cpp:65
const Geometry::Type geom
void GetBoundingBox(Vector &min, Vector &max, int ref=2)
Definition: mesh.cpp:95
RefinedGeometry * Refine(int Geom, int Times, int ETimes=1)
Definition: geom.cpp:802
int Size() const
Returns the size of the vector.
Definition: vector.hpp:113
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:614
STL namespace.
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:225
int material(Vector &x, Vector &xmin, Vector &xmax)
Definition: shaper.cpp:49
double Normlp(double p) const
Returns the l_p norm of the vector.
Definition: vector.cpp:727
int dim
Definition: ex3.cpp:47
void SetAttributes()
Definition: mesh.cpp:869
int Append(const T &el)
Append element to array, resize if necessary.
Definition: array.hpp:548
void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:3315
int MeshGenerator()
Get the mesh generator/type.
Definition: mesh.hpp:607
GeometryRefiner GlobGeometryRefiner
Definition: geom.cpp:1188
IntegrationRule RefPts
Definition: geom.hpp:210
int Dimension() const
Definition: mesh.hpp:641
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
int SpaceDimension() const
Definition: mesh.hpp:642
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:74
int GetElementBaseGeometry(int i=0) const
Definition: mesh.hpp:684
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:174
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
Definition: mesh.cpp:255
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:852
Vector data type.
Definition: vector.hpp:41
virtual void Transform(const IntegrationPoint &, Vector &)=0
void GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
Definition: mesh.cpp:6319
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:195
bool Good() const
Definition: optparser.hpp:120