libMesh
miscellaneous_ex10.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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 10 - Stitching meshes </h1>
21 // \author Dana Christen
22 // \date 2014
23 //
24 // This example shows how to generate a domain by stitching 8 cubic meshes
25 // together. Then a Poisson problem is solved on the stitched domain,
26 // and compared to the Poisson problem on a reference (unstitched) mesh
27 // to verify that the results match.
28 
29 
30 // C++ include files that we need
31 #include <iostream>
32 #include <algorithm>
33 #include <math.h>
34 #include <set>
35 #include <sstream>
36 #include <fstream>
37 
38 // libMesh includes
39 #include "libmesh/libmesh.h"
40 #include "libmesh/replicated_mesh.h"
41 #include "libmesh/mesh_generation.h"
42 #include "libmesh/linear_implicit_system.h"
43 #include "libmesh/equation_systems.h"
44 #include "libmesh/exact_solution.h"
45 #include "libmesh/exodusII_io.h"
46 #include "libmesh/fe.h"
47 #include "libmesh/quadrature_gauss.h"
48 #include "libmesh/dof_map.h"
49 #include "libmesh/sparse_matrix.h"
50 #include "libmesh/numeric_vector.h"
51 #include "libmesh/dense_matrix.h"
52 #include "libmesh/dense_vector.h"
53 #include "libmesh/elem.h"
54 #include "libmesh/dirichlet_boundaries.h"
55 #include "libmesh/zero_function.h"
56 #include "libmesh/libmesh_logging.h"
57 #include "libmesh/getpot.h"
58 #include "libmesh/error_vector.h"
59 #include "libmesh/kelly_error_estimator.h"
60 #include "libmesh/mesh_refinement.h"
61 
62 using namespace libMesh;
63 
64 bool compare_elements(const ReplicatedMesh & mesh1,
65  const ReplicatedMesh & mesh2);
67  const std::string & system_name);
69  EquationSystems &);
70 
71 int main (int argc, char ** argv)
72 {
73  START_LOG("Initialize and create cubes", "main");
74  LibMeshInit init (argc, argv);
75 
76  // This example requires a linear solver package.
77  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
78  "--enable-petsc, --enable-trilinos, or --enable-eigen");
79 
80  // Create a GetPot object to parse the command line
81  GetPot command_line (argc, argv);
82 
83  // Check for proper calling arguments.
84  if (argc < 3)
85  libmesh_error_msg("Usage:\n" << "\t " << argv[0] << " -n 15");
86 
87  // Brief message to the user regarding the program name
88  // and command line arguments.
89  else
90  {
91  libMesh::out << "Running " << argv[0];
92 
93  for (int i=1; i<argc; i++)
94  libMesh::out << " " << argv[i];
95 
96  libMesh::out << std::endl << std::endl;
97  }
98 
99  // This is 3D-only problem
100  const unsigned int dim = 3;
101 
102  // Skip higher-dimensional examples on a lower-dimensional libMesh build
103  libmesh_example_requires(dim <= LIBMESH_DIM, "3D support");
104 
105  // Read number of elements used in each cube from command line
106  int ps = 10;
107  if (command_line.search(1, "-n"))
108  ps = command_line.next(ps);
109 
110  // Generate eight meshes that will be stitched
111  ReplicatedMesh mesh (init.comm());
112  ReplicatedMesh mesh1(init.comm());
113  ReplicatedMesh mesh2(init.comm());
114  ReplicatedMesh mesh3(init.comm());
115  ReplicatedMesh mesh4(init.comm());
116  ReplicatedMesh mesh5(init.comm());
117  ReplicatedMesh mesh6(init.comm());
118  ReplicatedMesh mesh7(init.comm());
119  MeshTools::Generation::build_cube (mesh, ps, ps, ps, -1, 0, 0, 1, 0, 1, HEX8);
120  MeshTools::Generation::build_cube (mesh1, ps, ps, ps, 0, 1, 0, 1, 0, 1, HEX8);
121  MeshTools::Generation::build_cube (mesh2, ps, ps, ps, -1, 0, -1, 0, 0, 1, HEX8);
122  MeshTools::Generation::build_cube (mesh3, ps, ps, ps, 0, 1, -1, 0, 0, 1, HEX8);
123  MeshTools::Generation::build_cube (mesh4, ps, ps, ps, -1, 0, 0, 1, -1, 0, HEX8);
124  MeshTools::Generation::build_cube (mesh5, ps, ps, ps, 0, 1, 0, 1, -1, 0, HEX8);
125  MeshTools::Generation::build_cube (mesh6, ps, ps, ps, -1, 0, -1, 0, -1, 0, HEX8);
126  MeshTools::Generation::build_cube (mesh7, ps, ps, ps, 0, 1, -1, 0, -1, 0, HEX8);
127 
128  // Generate a single unstitched reference mesh
129  ReplicatedMesh nostitch_mesh(init.comm());
130  MeshTools::Generation::build_cube (nostitch_mesh, ps*2, ps*2, ps*2, -1, 1, -1, 1, -1, 1, HEX8);
131  STOP_LOG("Initialize and create cubes", "main");
132 
133  START_LOG("Stitching", "main");
134  // We stitch the meshes in a hierarchical way.
135  mesh.stitch_meshes(mesh1, 2, 4, TOLERANCE, true, true, false, false);
136  mesh2.stitch_meshes(mesh3, 2, 4, TOLERANCE, true, true, false, false);
137  mesh.stitch_meshes(mesh2, 1, 3, TOLERANCE, true, true, false, false);
138  mesh4.stitch_meshes(mesh5, 2, 4, TOLERANCE, true, true, false, false);
139  mesh6.stitch_meshes(mesh7, 2, 4, TOLERANCE, true, true, false, false);
140  mesh4.stitch_meshes(mesh6, 1, 3, TOLERANCE, true, true, false, false);
141  mesh.stitch_meshes(mesh4, 0, 5, TOLERANCE, true, true, false, false);
142  STOP_LOG("Stitching", "main");
143 
144  START_LOG("Initialize and solve systems", "main");
145  EquationSystems equation_systems_stitch (mesh);
146  assemble_and_solve(mesh, equation_systems_stitch);
147 
148  EquationSystems equation_systems_nostitch (nostitch_mesh);
149  assemble_and_solve(nostitch_mesh, equation_systems_nostitch);
150  STOP_LOG("Initialize and solve systems", "main");
151 
152  START_LOG("Result comparison", "main");
153  ExactSolution comparison(equation_systems_stitch);
154  comparison.attach_reference_solution(&equation_systems_nostitch);
155  comparison.compute_error("Poisson", "u");
156  Real error = comparison.l2_error("Poisson", "u");
157  libmesh_assert_less(error, TOLERANCE*sqrt(TOLERANCE));
158  libMesh::out << "L2 error between stitched and non-stitched cases: " << error << std::endl;
159  STOP_LOG("Result comparison", "main");
160 
161  START_LOG("Output", "main");
162 #ifdef LIBMESH_HAVE_EXODUS_API
163  ExodusII_IO(mesh).write_equation_systems("solution_stitch.exo",
164  equation_systems_stitch);
165  ExodusII_IO(nostitch_mesh).write_equation_systems("solution_nostitch.exo",
166  equation_systems_nostitch);
167 #endif // #ifdef LIBMESH_HAVE_EXODUS_API
168  STOP_LOG("Output", "main");
169 
170  return 0;
171 }
172 
174  EquationSystems & equation_systems)
175 {
176  mesh.print_info();
177 
178  LinearImplicitSystem & system =
179  equation_systems.add_system<LinearImplicitSystem> ("Poisson");
180 
181  unsigned int u_var = system.add_variable("u", FIRST, LAGRANGE);
182 
184 
185  // the cube has boundaries IDs 0, 1, 2, 3, 4 and 5
186  std::set<boundary_id_type> boundary_ids;
187  for (int j = 0; j<6; ++j)
188  boundary_ids.insert(j);
189 
190  // Create a vector storing the variable numbers which the BC applies to
191  std::vector<unsigned int> variables(1);
192  variables[0] = u_var;
193 
194  ZeroFunction<> zf;
195 
196  // Most DirichletBoundary users will want to supply a "locally
197  // indexed" functor
198  DirichletBoundary dirichlet_bc(boundary_ids, variables, zf,
200  system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);
201 
202  equation_systems.init();
203  equation_systems.print_info();
204 
205 #ifdef LIBMESH_ENABLE_AMR
206  MeshRefinement mesh_refinement(mesh);
207 
208  mesh_refinement.refine_fraction() = 0.7;
209  mesh_refinement.coarsen_fraction() = 0.3;
210  mesh_refinement.max_h_level() = 5;
211 
212  const unsigned int max_r_steps = 2;
213 
214  for (unsigned int r_step=0; r_step<=max_r_steps; r_step++)
215  {
216  system.solve();
217  if (r_step != max_r_steps)
218  {
219  ErrorVector error;
220  KellyErrorEstimator error_estimator;
221 
222  error_estimator.estimate_error(system, error);
223 
224  libMesh::out << "Error estimate\nl2 norm = "
225  << error.l2_norm()
226  << "\nmaximum = "
227  << error.maximum()
228  << std::endl;
229 
230  mesh_refinement.flag_elements_by_error_fraction (error);
231 
232  mesh_refinement.refine_and_coarsen_elements();
233 
234  equation_systems.reinit();
235  }
236  }
237 #else
238  system.solve();
239 #endif
240 }
241 
243  const std::string & libmesh_dbg_var(system_name))
244 {
245  libmesh_assert_equal_to (system_name, "Poisson");
246 
247  const MeshBase & mesh = es.get_mesh();
248  const unsigned int dim = mesh.mesh_dimension();
249  LinearImplicitSystem & system = es.get_system<LinearImplicitSystem>("Poisson");
250 
251  const DofMap & dof_map = system.get_dof_map();
252 
253  FEType fe_type = dof_map.variable_type(0);
254  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
255  QGauss qrule (dim, FIFTH);
256  fe->attach_quadrature_rule (&qrule);
257  UniquePtr<FEBase> fe_face (FEBase::build(dim, fe_type));
258  QGauss qface(dim-1, FIFTH);
259  fe_face->attach_quadrature_rule (&qface);
260 
261  const std::vector<Real> & JxW = fe->get_JxW();
262  const std::vector<std::vector<Real>> & phi = fe->get_phi();
263  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
264 
267 
268  std::vector<dof_id_type> dof_indices;
269 
270  for (const auto & elem : mesh.active_local_element_ptr_range())
271  {
272  dof_map.dof_indices (elem, dof_indices);
273 
274  fe->reinit (elem);
275 
276  Ke.resize (dof_indices.size(),
277  dof_indices.size());
278 
279  Fe.resize (dof_indices.size());
280 
281  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
282  {
283  for (std::size_t i=0; i<phi.size(); i++)
284  {
285  Fe(i) += JxW[qp]*phi[i][qp];
286  for (std::size_t j=0; j<phi.size(); j++)
287  Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
288  }
289  }
290 
291  dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
292 
293  system.matrix->add_matrix (Ke, dof_indices);
294  system.rhs->add_vector (Fe, dof_indices);
295  }
296 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
This class handles the computation of the L2 and/or H1 error for the Systems in the EquationSystems o...
This is the EquationSystems class.
The ReplicatedMesh class is derived from the MeshBase class, and is used to store identical copies of...
bool compare_elements(const ReplicatedMesh &mesh1, const ReplicatedMesh &mesh2)
ConstFunction that simply returns 0.
Definition: zero_function.h:35
unsigned int dim
The ErrorVector is a specialization of the StatisticsVector for error data computed on a finite eleme...
Definition: error_vector.h:50
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
Definition: exodusII_io.h:52
This class provides a specific system class.
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:350
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]...
virtual Real l2_norm() const
Definition: statistics.C:36
virtual T maximum() const
Definition: statistics.C:61
unsigned int add_variable(const std::string &var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=libmesh_nullptr)
Adds the variable var to the list of variables for this system.
Definition: system.C:1101
Real & coarsen_fraction()
The coarsen_fraction sets either a desired target or a desired maximum number of elements to flag for...
MeshBase & mesh
NumericVector< Number > * rhs
The system matrix.
This class allows one to associate Dirichlet boundary values with a given set of mesh boundary ids an...
Real & refine_fraction()
The refine_fraction sets either a desired target or a desired maximum number of elements to flag for ...
void assemble_and_solve(MeshBase &, EquationSystems &)
static const Real TOLERANCE
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
Definition: libmesh.h:62
The libMesh namespace provides an interface to certain functionality in the library.
unsigned int & max_h_level()
The max_h_level is the greatest refinement level an element should reach.
void compute_error(const std::string &sys_name, const std::string &unknown_name)
Computes and stores the error in the solution value e = u-u_h, the gradient grad(e) = grad(u) - grad(...
This is the MeshBase class.
Definition: mesh_base.h:68
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
SolverPackage default_solver_package()
Definition: libmesh.C:995
virtual SimpleRange< element_iterator > active_local_element_ptr_range()=0
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
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.
const DofMap & get_dof_map() const
Definition: system.h:2030
This is the MeshRefinement class.
void assemble_poisson(EquationSystems &es, const std::string &system_name)
void print_info(std::ostream &os=libMesh::out) const
Prints information about the equation systems, by default to libMesh::out.
void attach_assemble_function(void fptr(EquationSystems &es, const std::string &name))
Register a user function to use in assembling the system matrix and RHS.
Definition: system.C:1795
virtual void reinit()
Reinitialize all the systems.
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
virtual System & add_system(const std::string &system_type, const std::string &name)
Add the system of type system_type named name to the systems array.
int main(int argc, char **argv)
const Parallel::Communicator & comm() const
Definition: libmesh.h:81
PetscErrorCode Vec Mat libmesh_dbg_var(j)
virtual void write_equation_systems(const std::string &, const EquationSystems &, const std::set< std::string > *system_names=libmesh_nullptr)
This method implements writing a mesh with data to a specified file where the data is taken from the ...
Definition: mesh_output.C:31
This class implements the Kelly error indicator which is based on the flux jumps between elements...
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
OStreamProxy out
SparseMatrix< Number > * matrix
The system matrix.
unsigned int mesh_dimension() const
Definition: mesh_base.C:148
bool refine_and_coarsen_elements()
Refines and coarsens user-requested elements.
virtual void estimate_error(const System &system, ErrorVector &error_per_cell, const NumericVector< Number > *solution_vector=libmesh_nullptr, bool estimate_parent_error=false) libmesh_override
This function uses the derived class&#39;s jump error estimate formula to estimate the error on each cell...
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:776
This class implements specific orders of Gauss quadrature.
const MeshBase & get_mesh() const
void flag_elements_by_error_fraction(const ErrorVector &error_per_cell, const Real refine_fraction=0.3, const Real coarsen_fraction=0.0, const unsigned int max_level=libMesh::invalid_uint)
Flags elements for coarsening and refinement based on the computed error passed in error_per_cell...
void add_dirichlet_boundary(const DirichletBoundary &dirichlet_boundary)
Adds a copy of the specified Dirichlet boundary to the system.
const T_sys & get_system(const std::string &name) const
virtual void init()
Initialize all the systems.
void attach_reference_solution(const EquationSystems *es_fine)
Attach function similar to system.h which allows the user to attach a second EquationSystems object w...
Real l2_error(const std::string &sys_name, const std::string &unknown_name)
virtual void solve() libmesh_override
Assembles & solves the linear system A*x=b.
void print_info(std::ostream &os=libMesh::out) const
Prints relevant information about the mesh.
Definition: mesh_base.C:448
unsigned int n_points() const
Definition: quadrature.h:113
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.
static UniquePtr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.