libMesh
miscellaneous_ex1.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // <h1>Miscellaneous Example 1 - Infinite Elements for the Wave Equation</h1>
21 // \author Daniel Dreyer
22 // \date 2003
23 //
24 // This is the sixth example program. It builds on
25 // the previous examples, and introduces the Infinite
26 // Element class. Note that the library must be compiled
27 // with Infinite Elements enabled. Otherwise, this
28 // example will abort.
29 // This example intends to demonstrate the similarities
30 // between the FE and the InfFE classes in libMesh.
31 // The matrices are assembled according to the wave equation.
32 // However, for practical applications a time integration
33 // scheme (as introduced in subsequent examples) should be
34 // used.
35 
36 // C++ include files that we need
37 #include <iostream>
38 #include <algorithm>
39 #include <math.h>
40 
41 // Basic include files needed for the mesh functionality.
42 #include "libmesh/exodusII_io.h"
43 #include "libmesh/libmesh.h"
44 #include "libmesh/mesh.h"
45 #include "libmesh/mesh_generation.h"
46 #include "libmesh/linear_implicit_system.h"
47 #include "libmesh/equation_systems.h"
48 #include "libmesh/enum_xdr_mode.h"
49 #include "libmesh/getpot.h"
50 #include "libmesh/mesh_refinement.h"
51 
52 // Define the Finite and Infinite Element object.
53 #include "libmesh/fe.h"
54 #include "libmesh/inf_fe.h"
55 #include "libmesh/inf_elem_builder.h"
56 
57 // Define Gauss quadrature rules.
58 #include "libmesh/quadrature_gauss.h"
59 
60 // Define useful datatypes for finite element
61 // matrix and vector components.
62 #include "libmesh/sparse_matrix.h"
63 #include "libmesh/numeric_vector.h"
64 #include "libmesh/dense_matrix.h"
65 #include "libmesh/dense_vector.h"
66 
67 // Define the DofMap, which handles degree of freedom
68 // indexing.
69 #include "libmesh/dof_map.h"
70 
71 // The definition of a vertex associated with a Mesh.
72 #include "libmesh/node.h"
73 
74 // The definition of a geometric element
75 #include "libmesh/elem.h"
76 
77 // Bring in everything from the libMesh namespace
78 using namespace libMesh;
79 
80 // Function prototype. This is similar to the Poisson
81 // assemble function of example 4.
83  const std::string & system_name);
84 
85 // Begin the main program.
86 int main (int argc, char ** argv)
87 {
88  // Initialize libMesh, like in example 2.
89  LibMeshInit init (argc, argv);
90 
91  // This example requires Infinite Elements
92 #ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
93  libmesh_example_requires(false, "--enable-ifem");
94 #else
95 
96  // Skip this 3D example if libMesh was compiled as 1D/2D-only.
97  libmesh_example_requires(3 <= LIBMESH_DIM, "3D support");
98 
99  // Create a serialized mesh, distributed across the default MPI
100  // communicator.
101  Mesh mesh(init.comm());
102 
103  // Get command line arguments for mesh size
104  GetPot input(argc, argv);
105 
106  const unsigned int nx = input("nx", 4),
107  ny = input("ny", 4),
108  nz = input("nz", 4);
109 
110  // Use the internal mesh generator to create elements
111  // on the square [-1,1]^3, of type Hex8.
113  nx, ny, nz,
114  -1., 1.,
115  -1., 1.,
116  -1., 1.,
117  HEX8);
118 
119  // Print information about the mesh to the screen.
120  mesh.print_info();
121 
122  // Write the mesh before the infinite elements are added
123 #ifdef LIBMESH_HAVE_EXODUS_API
124  ExodusII_IO(mesh).write ("orig_mesh.e");
125 #endif
126 
127  // Normally, when a mesh is imported or created in
128  // libMesh, only conventional elements exist. The infinite
129  // elements used here, however, require prescribed
130  // nodal locations (with specified distances from an imaginary
131  // origin) and configurations that a conventional mesh creator
132  // in general does not offer. Therefore, an efficient method
133  // for building infinite elements is offered. It can account
134  // for symmetry planes and creates infinite elements in a fully
135  // automatic way.
136  //
137  // Right now, the simplified interface is used, automatically
138  // determining the origin. Check MeshBase for a generalized
139  // method that can even return the element faces of interior
140  // vibrating surfaces. The bool determines whether to be
141  // verbose.
142  InfElemBuilder builder(mesh);
143  builder.build_inf_elem(true);
144 
145  // Reassign subdomain_id() of all infinite elements.
146  // Otherwise, the exodus-api will fail on the mesh.
147  for (auto & elem : mesh.element_ptr_range())
148  if (elem->infinite())
149  elem->subdomain_id() = 1;
150 
151  // Print information about the mesh to the screen.
152  mesh.print_info();
153 
154  // Write the mesh with the infinite elements added.
155  // Compare this to the original mesh.
156 #ifdef LIBMESH_HAVE_EXODUS_API
157  ExodusII_IO(mesh).write ("ifems_added.e");
158 #endif
159 
160  // After building infinite elements, we have to let
161  // the elements find their neighbors again.
163 
164  // Create an equation systems object, where ThinSystem
165  // offers only the crucial functionality for solving a
166  // system. Use ThinSystem when you want the sleekest
167  // system possible.
168  EquationSystems equation_systems (mesh);
169 
170  // Declare the system and its variables.
171  // Create a system named "Wave". This can
172  // be a simple, steady system
173  equation_systems.add_system<LinearImplicitSystem> ("Wave");
174 
175  // Create an FEType describing the approximation
176  // characteristics of the InfFE object. Note that
177  // the constructor automatically defaults to some
178  // sensible values. But use FIRST order
179  // approximation.
180  FEType fe_type(FIRST);
181 
182  // Add the variable "p" to "Wave". Note that there exist
183  // various approaches in adding variables. In example 3,
184  // add_variable took the order of approximation and used
185  // default values for the FEFamily, while here the FEType
186  // is used.
187  equation_systems.get_system("Wave").add_variable("p", fe_type);
188 
189  // Give the system a pointer to the matrix assembly
190  // function.
191  equation_systems.get_system("Wave").attach_assemble_function (assemble_wave);
192 
193  // Set the speed of sound and fluid density
194  // as EquationSystems parameter,
195  // so that assemble_wave() can access it.
196  equation_systems.parameters.set<Real>("speed") = 1.;
197  equation_systems.parameters.set<Real>("fluid density") = 1.;
198 
199  // Initialize the data structures for the equation system.
200  equation_systems.init();
201 
202 #ifdef LIBMESH_ENABLE_AMR
203  // Do uniform refinement if requested
204  const unsigned int nr = input("nr", 0);
205  if (nr)
206  {
207  MeshRefinement mesh_refinement(mesh);
208  mesh_refinement.uniformly_refine(nr);
209  equation_systems.reinit();
210  equation_systems.print_info();
211  }
212 #endif
213 
214  // Print and solve the refined sysem
215  equation_systems.get_system("Wave").solve();
216 
217  libMesh::out << "Wave system solved" << std::endl;
218 
219  // Write the whole EquationSystems object to file.
220  // For infinite elements, the concept of nodal_soln()
221  // is not applicable. Therefore, writing the mesh in
222  // some format @e always gives all-zero results at
223  // the nodes of the infinite elements. Instead,
224  // use the FEInterface::compute_data() methods to
225  // determine physically correct results within an
226  // infinite element.
227  equation_systems.write ("eqn_sys.dat", WRITE);
228 
229  libMesh::out << "eqn_sys.dat written" << std::endl;
230 
231  // All done.
232  return 0;
233 
234 #endif // else part of ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
235 }
236 
237 // This function assembles the system matrix and right-hand-side
238 // for the discrete form of our wave equation.
240  const std::string & libmesh_dbg_var(system_name))
241 {
242  // It is a good idea to make sure we are assembling
243  // the proper system.
244  libmesh_assert_equal_to (system_name, "Wave");
245 
246  // Avoid unused variable warnings when compiling without infinite
247  // elements enabled.
248  libmesh_ignore(es);
249 
250 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
251 
252  // Get a constant reference to the mesh object.
253  const MeshBase & mesh = es.get_mesh();
254 
255  // Get a reference to the system we are solving.
257 
258  // A reference to the DofMap object for this system. The DofMap
259  // object handles the index translation from node and element numbers
260  // to degree of freedom numbers.
261  const DofMap & dof_map = system.get_dof_map();
262 
263  // The dimension that we are running.
264  const unsigned int dim = mesh.mesh_dimension();
265 
266  // Copy the speed of sound to a local variable.
267  const Real speed = es.parameters.get<Real>("speed");
268 
269  // Get a constant reference to the Finite Element type
270  // for the first (and only) variable in the system.
271  const FEType & fe_type = dof_map.variable_type(0);
272 
273  // Build a Finite Element object of the specified type. Since the
274  // FEBase::build() member dynamically creates memory we will
275  // store the object as a std::unique_ptr<FEBase>.
276  std::unique_ptr<FEBase> fe (FEBase::build(dim, fe_type));
277 
278  // Do the same for an infinite element.
279  std::unique_ptr<FEBase> inf_fe (FEBase::build_InfFE(dim, fe_type));
280 
281  // A 2nd order Gauss quadrature rule for numerical integration.
282  QGauss qrule (dim, SECOND);
283 
284  // Tell the finite element object to use our quadrature rule.
285  fe->attach_quadrature_rule (&qrule);
286 
287  // Due to its internal structure, the infinite element handles
288  // quadrature rules differently. It takes the quadrature
289  // rule which has been initialized for the FE object, but
290  // creates suitable quadrature rules by @e itself. The user
291  // need not worry about this.
292  inf_fe->attach_quadrature_rule (&qrule);
293 
294  // Define data structures to contain the element matrix
295  // and right-hand-side vector contribution. Following
296  // basic finite element terminology we will denote these
297  // "Ke", "Ce", "Me", and "Fe" for the stiffness, damping
298  // and mass matrices, and the load vector. Note that in
299  // Acoustics, these descriptors though do @e not match the
300  // true physical meaning of the projectors. The final
301  // overall system, however, resembles the conventional
302  // notation again.
307 
308  // This vector will hold the degree of freedom indices for
309  // the element. These define where in the global system
310  // the element degrees of freedom get mapped.
311  std::vector<dof_id_type> dof_indices;
312 
313  // The global system matrix
314  SparseMatrix<Number> & matrix = system.get_system_matrix();
315 
316  // Now we will loop over all the elements in the mesh.
317  // We will compute the element matrix and right-hand-side
318  // contribution.
319  for (const auto & elem : mesh.active_local_element_ptr_range())
320  {
321  // Get the degree of freedom indices for the
322  // current element. These define where in the global
323  // matrix and right-hand-side this element will
324  // contribute to.
325  dof_map.dof_indices (elem, dof_indices);
326 
327  const unsigned int n_dofs =
328  cast_int<unsigned int>(dof_indices.size());
329 
330  // The mesh contains both finite and infinite elements. These
331  // elements are handled through different classes, namely
332  // FE and InfFE, respectively. However, since both
333  // are derived from FEBase, they share the same interface,
334  // and overall burden of coding is @e greatly reduced through
335  // using a pointer, which is adjusted appropriately to the
336  // current element type.
337  FEBase * cfe = nullptr;
338 
339  // This here is almost the only place where we need to
340  // distinguish between finite and infinite elements.
341  // For faster computation, however, different approaches
342  // may be feasible.
343  //
344  // Up to now, we do not know what kind of element we
345  // have. Aske the element of what type it is:
346  if (elem->infinite())
347  {
348  // We have an infinite element. Let cfe point
349  // to our InfFE object. This is handled through
350  // a std::unique_ptr. Through the std::unique_ptr::get() we "borrow"
351  // the pointer, while the std::unique_ptr inf_fe is
352  // still in charge of memory management.
353  cfe = inf_fe.get();
354  }
355  else
356  {
357  // This is a conventional finite element. Let fe handle it.
358  cfe = fe.get();
359 
360  // Boundary conditions.
361  // Here we just zero the rhs-vector. For natural boundary
362  // conditions check e.g. previous examples.
363  {
364  // Zero the RHS for this element.
365  Fe.resize (n_dofs);
366 
367  system.rhs->add_vector (Fe, dof_indices);
368  } // end boundary condition section
369  } // else if (elem->infinite())
370 
371  // This is slightly different from the Poisson solver:
372  // Since the finite element object may change, we have to
373  // initialize the constant references to the data fields
374  // each time again, when a new element is processed.
375  //
376  // Due to infinite extent of the element, the 'pure' Jacobian is divergent
377  // and we must use a re-scaled one. As re-scaling, a 'decay'-function
378  // is used (in 3D this is \f$ r^{-2} \f$).
379  //
380  // To account for this extra weight, \p phi, \p dphi and \p weight are
381  // re-scaled as well:
382  // * J --> J x decay^2
383  // * phi --> phi/decay x r
384  // * dphi --> dphi/decay x r
385  // * weight --> weight / r^2
386  // With this, the product of Jacobian, test and trial functions (and their derivatives)
387  // can be used as normal since the extra weights cancel.
388  //
389  // The element Jacobian * quadrature weight at each integration point.
390  const std::vector<Real> & JxW = cfe->get_JxWxdecay_sq();
391 
392  // The element shape functions evaluated at the quadrature points.
393  const std::vector<std::vector<Real>> & phi = cfe->get_phi_over_decayxR();
394  //const std::vector<std::vector<Real>> & phi = cfe->get_phi();
395 
396  // The element shape function gradients evaluated at the quadrature
397  // points, divided by weight x rad
398  const std::vector<std::vector<RealGradient>> & dphi = cfe->get_dphi_over_decayxR();
399  //const std::vector<std::vector<RealGradient>> & dphi = cfe->get_dphi();
400 
401  // The infinite elements need more data fields than conventional FE.
402  // These are the gradients of the phase term dphase, an additional
403  // radial weight for the test functions Sobolev_weight, and its
404  // gradient.
405  //
406  // Note that these data fields are also initialized appropriately by
407  // the FE method, so that the weak form (below) is valid for @e both
408  // finite and infinite elements.
409  const std::vector<RealGradient> & dphase = cfe->get_dphase();
410  const std::vector<Real> & weight = cfe->get_Sobolev_weightxR_sq();
411  const std::vector<RealGradient> & dweight = cfe->get_Sobolev_dweightxR_sq();
412 
413  // Now this is all independent of whether we use an FE
414  // or an InfFE. Nice, hm? ;-)
415  //
416  // Compute the element-specific data, as described
417  // in previous examples.
418  cfe->reinit (elem);
419 
420  // Zero the element matrices. Boundary conditions were already
421  // processed in the FE-only section, see above.
422  Ke.resize (n_dofs, n_dofs);
423  Ce.resize (n_dofs, n_dofs);
424  Me.resize (n_dofs, n_dofs);
425 
426  // The total number of quadrature points for infinite elements
427  // @e has to be determined in a different way, compared to
428  // conventional finite elements. This type of access is also
429  // valid for finite elements, so this can safely be used
430  // anytime, instead of asking the quadrature rule, as
431  // seen in previous examples.
432  unsigned int max_qp = cfe->n_quadrature_points();
433 
434  // Loop over the quadrature points.
435  for (unsigned int qp=0; qp<max_qp; qp++)
436  {
437  // Similar to the modified access to the number of quadrature
438  // points, the number of shape functions may also be obtained
439  // in a different manner. This offers the great advantage
440  // of being valid for both finite and infinite elements.
441  const unsigned int n_sf =
442  FEInterface::n_dofs(cfe->get_fe_type(), elem);
443 
444  // Now we will build the element matrices. Since the infinite
445  // elements are based on a Petrov-Galerkin scheme, the
446  // resulting system matrices are non-symmetric. The additional
447  // weight, described before, is part of the trial space.
448  //
449  // For the finite elements, though, these matrices are symmetric
450  // just as we know them, since the additional fields dphase,
451  // weight, and dweight are initialized appropriately.
452  //
453  // test functions: weight[qp]*phi[i][qp]
454  // trial functions: phi[j][qp]
455  // phase term: phase[qp]
456  //
457  // derivatives are similar, but note that these are of type
458  // Point, not of type Real.
459  for (unsigned int i=0; i<n_sf; i++)
460  for (unsigned int j=0; j<n_sf; j++)
461  {
462  // (ndt*Ht + nHt*d) * nH
463  Ke(i,j) +=
464  (
465  (dweight[qp] * phi[i][qp] // Point * Real = Point
466  + // +
467  dphi[i][qp] * weight[qp] // Point * Real = Point
468  ) * dphi[j][qp]
469  ) * JxW[qp];
470 
471  // (d*Ht*nmut*nH - ndt*nmu*Ht*H - d*nHt*nmu*H)
472  Ce(i,j) +=
473  (
474  (dphase[qp] * dphi[j][qp]) // (Point * Point) = Real
475  * weight[qp] * phi[i][qp] // * Real * Real = Real
476  - // -
477  (dweight[qp] * dphase[qp]) // (Point * Point) = Real
478  * phi[i][qp] * phi[j][qp] // * Real * Real = Real
479  - // -
480  (dphi[i][qp] * dphase[qp]) // (Point * Point) = Real
481  * weight[qp] * phi[j][qp] // * Real * Real = Real
482  ) * JxW[qp];
483 
484  // (d*Ht*H * (1 - nmut*nmu))
485  Me(i,j) +=
486  (
487  (1. - (dphase[qp] * dphase[qp])) // (Real - (Point * Point)) = Real
488  * phi[i][qp] * phi[j][qp] * weight[qp] // * Real * Real * Real = Real
489  ) * JxW[qp];
490 
491  } // end of the matrix summation loop
492  } // end of quadrature point loop
493 
494  // The element matrices are now built for this element.
495  // Collect them in Ke, and then add them to the global matrix.
496  // The SparseMatrix::add_matrix() member does this for us.
497  Ke.add(1./speed , Ce);
498  Ke.add(1./(speed*speed), Me);
499 
500  // If this assembly program were to be used on an adaptive mesh,
501  // we would have to apply any hanging node constraint equations
502  dof_map.constrain_element_matrix(Ke, dof_indices);
503 
504  matrix.add_matrix (Ke, dof_indices);
505  } // end of element loop
506 
507  // Note that we have not applied any boundary conditions so far.
508  // Here we apply a unit load at the node located at (0,0,0).
509  for (const auto & node : mesh.local_node_ptr_range())
510  if (std::abs((*node)(0)) < TOLERANCE &&
511  std::abs((*node)(1)) < TOLERANCE &&
512  std::abs((*node)(2)) < TOLERANCE)
513  {
514  // The global number of the respective degree of freedom.
515  unsigned int dn = node->dof_number(0,0,0);
516 
517  system.rhs->add (dn, 1.);
518  }
519 
520 #else
521 
522  // dummy assert
523  libmesh_assert_not_equal_to (es.get_mesh().mesh_dimension(), 1);
524 
525 #endif //ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
526 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:196
This is the EquationSystems class.
int main(int argc, char **argv)
static unsigned int n_dofs(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface.C:355
void write(std::string_view name, const XdrMODE, const unsigned int write_flags=(WRITE_DATA), bool partition_agnostic=true) const
Write the systems to disk using the XDR data format.
static constexpr Real TOLERANCE
unsigned int dim
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
Definition: exodusII_io.h:52
Manages consistently variables, degrees of freedom, coefficient vectors, matrices and linear solvers ...
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:396
virtual void add_vector(const T *v, const std::vector< numeric_index_type > &dof_indices)
Computes , where v is a pointer and each dof_indices[i] specifies where to add value v[i]...
void print_info(std::ostream &os=libMesh::out) const
Prints information about the equation systems, by default to libMesh::out.
MeshBase & mesh
NumericVector< Number > * rhs
The system matrix.
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
Definition: libmesh.h:90
The libMesh namespace provides an interface to certain functionality in the library.
const T_sys & get_system(std::string_view name) const
This is the MeshBase class.
Definition: mesh_base.h:75
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:179
void libmesh_ignore(const Args &...)
virtual void add_matrix(const DenseMatrix< T > &dm, const std::vector< numeric_index_type > &rows, const std::vector< numeric_index_type > &cols)=0
Add the full matrix dm to the SparseMatrix.
virtual void find_neighbors(const bool reset_remote_elements=false, const bool reset_current_list=true)=0
Locate element face (edge in 2D) neighbors.
const T & get(std::string_view) const
Definition: parameters.h:426
Implements (adaptive) mesh refinement algorithms for a MeshBase.
dof_id_type weight(const MeshBase &mesh, const processor_id_type pid)
Definition: mesh_tools.C:437
void print_info(std::ostream &os=libMesh::out, const unsigned int verbosity=0, const bool global=true) const
Prints relevant information about the mesh.
Definition: mesh_base.C:1562
virtual void reinit()
Handle any mesh changes and reinitialize all the systems on the updated mesh.
boostcopy::enable_if_c< ScalarTraits< T2 >::value, void >::type add(const T2 factor, const DenseMatrix< T3 > &mat)
Adds factor times mat to this matrix.
void init(triangulateio &t)
Initializes the fields of t to nullptr/0 as necessary.
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
This class is used to build infinite elements on top of an existing mesh.
void assemble_wave(EquationSystems &es, const std::string &system_name)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const Point build_inf_elem(const bool be_verbose=false)
Build infinite elements atop a volume-based mesh, determine origin automatically. ...
virtual void write(const std::string &fname) override
This method implements writing a mesh to a specified file.
Definition: exodusII_io.C:2180
T & set(const std::string &)
Definition: parameters.h:469
OStreamProxy out
const MeshBase & get_mesh() const
void resize(const unsigned int new_m, const unsigned int new_n)
Resizes the matrix to the specified size and calls zero().
Definition: dense_matrix.h:895
This class implements specific orders of Gauss quadrature.
unsigned int mesh_dimension() const
Definition: mesh_base.C:372
Parameters parameters
Data structure holding arbitrary parameters.
static std::unique_ptr< FEGenericBase > build_InfFE(const unsigned int dim, const FEType &type)
Builds a specific infinite element type.
virtual void init()
Initialize all the systems.
virtual void add(const numeric_index_type i, const T value)=0
Adds value to the vector entry specified by i.
virtual System & add_system(std::string_view system_type, std::string_view name)
Add the system of type system_type named name to the systems array.
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
const DofMap & get_dof_map() const
Definition: system.h:2374
const SparseMatrix< Number > & get_system_matrix() const
void build_cube(UnstructuredMesh &mesh, const unsigned int nx=0, const unsigned int ny=0, const unsigned int nz=0, const Real xmin=0., const Real xmax=1., const Real ymin=0., const Real ymax=1., const Real zmin=0., const Real zmax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
Builds a (elements) cube.
This class forms the foundation from which generic finite elements may be derived.
void uniformly_refine(unsigned int n=1)
Uniformly refines the mesh n times.