libMesh
Functions
miscellaneous_ex10.C File Reference

Go to the source code of this file.

Functions

bool compare_elements (const ReplicatedMesh &mesh1, const ReplicatedMesh &mesh2)
 
void assemble_poisson (EquationSystems &es, const std::string &system_name)
 
void assemble_and_solve (MeshBase &, EquationSystems &)
 
int main (int argc, char **argv)
 
void assemble_poisson (EquationSystems &es, const std::string &libmesh_dbg_var(system_name))
 

Function Documentation

void assemble_and_solve ( MeshBase mesh,
EquationSystems equation_systems 
)

Definition at line 173 of file miscellaneous_ex10.C.

References libMesh::DofMap::add_dirichlet_boundary(), libMesh::EquationSystems::add_system(), libMesh::System::add_variable(), assemble_poisson(), libMesh::System::attach_assemble_function(), libMesh::MeshRefinement::coarsen_fraction(), libMesh::JumpErrorEstimator::estimate_error(), libMesh::FIRST, libMesh::MeshRefinement::flag_elements_by_error_fraction(), libMesh::System::get_dof_map(), libMesh::EquationSystems::init(), libMesh::StatisticsVector< T >::l2_norm(), libMesh::LAGRANGE, libMesh::LOCAL_VARIABLE_ORDER, libMesh::MeshRefinement::max_h_level(), libMesh::StatisticsVector< T >::maximum(), libMesh::out, libMesh::EquationSystems::print_info(), libMesh::MeshBase::print_info(), libMesh::MeshRefinement::refine_and_coarsen_elements(), libMesh::MeshRefinement::refine_fraction(), libMesh::EquationSystems::reinit(), and libMesh::LinearImplicitSystem::solve().

Referenced by main().

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 }
ConstFunction that simply returns 0.
Definition: zero_function.h:35
The ErrorVector is a specialization of the StatisticsVector for error data computed on a finite eleme...
Definition: error_vector.h:50
This class provides a specific system class.
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
This class allows one to associate Dirichlet boundary values with a given set of mesh boundary ids an...
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.
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.
This class implements the Kelly error indicator which is based on the flux jumps between elements...
OStreamProxy out
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 add_dirichlet_boundary(const DirichletBoundary &dirichlet_boundary)
Adds a copy of the specified Dirichlet boundary to the system.
virtual void init()
Initialize all the systems.
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
void assemble_poisson ( EquationSystems es,
const std::string &  system_name 
)

Referenced by assemble_and_solve().

void assemble_poisson ( EquationSystems es,
const std::string &  libmesh_dbg_varsystem_name 
)

Definition at line 242 of file miscellaneous_ex10.C.

References libMesh::MeshBase::active_local_element_ptr_range(), libMesh::SparseMatrix< T >::add_matrix(), libMesh::NumericVector< T >::add_vector(), libMesh::FEGenericBase< OutputType >::build(), dim, libMesh::FIFTH, libMesh::System::get_dof_map(), libMesh::EquationSystems::get_mesh(), libMesh::EquationSystems::get_system(), libMesh::ImplicitSystem::matrix, mesh, libMesh::MeshBase::mesh_dimension(), libMesh::QBase::n_points(), libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), and libMesh::ExplicitSystem::rhs.

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
unsigned int dim
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]...
MeshBase & mesh
NumericVector< Number > * rhs
The system matrix.
This is the MeshBase class.
Definition: mesh_base.h:68
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
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
SparseMatrix< Number > * matrix
The system matrix.
unsigned int mesh_dimension() const
Definition: mesh_base.C:148
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
const T_sys & get_system(const std::string &name) const
bool compare_elements ( const ReplicatedMesh mesh1,
const ReplicatedMesh mesh2 
)
int main ( int  argc,
char **  argv 
)

Definition at line 71 of file miscellaneous_ex10.C.

References assemble_and_solve(), libMesh::ExactSolution::attach_reference_solution(), libMesh::MeshTools::Generation::build_cube(), libMesh::LibMeshInit::comm(), libMesh::ExactSolution::compute_error(), libMesh::default_solver_package(), dim, libMesh::HEX8, libMesh::TriangleWrapper::init(), libMesh::INVALID_SOLVER_PACKAGE, libMesh::ExactSolution::l2_error(), mesh, libMesh::out, libMesh::Real, libMesh::TOLERANCE, and libMesh::MeshOutput< MT >::write_equation_systems().

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 }
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...
unsigned int dim
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
Definition: exodusII_io.h:52
MeshBase & mesh
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
SolverPackage default_solver_package()
Definition: libmesh.C:995
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
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
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
OStreamProxy out
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.