libMesh
Functions
vector_fe_ex4.C File Reference

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Function Documentation

◆ main()

int main ( int  argc,
char **  argv 
)

Definition at line 48 of file vector_fe_ex4.C.

References libMesh::DiffSolver::absolute_residual_tolerance, libMesh::EquationSystems::add_system(), libMesh::ExactSolution::attach_exact_derivs(), libMesh::ExactSolution::attach_exact_values(), libMesh::MeshTools::Generation::build_cube(), libMesh::command_line_value(), libMesh::ExactSolution::compute_error(), libMesh::default_solver_package(), libMesh::ExactSolution::error_norm(), libMesh::ExactSolution::extra_quadrature_order(), libMesh::ExactSolution::hcurl_error(), libMesh::HCURL_SEMINORM, libMesh::TriangleWrapper::init(), libMesh::EquationSystems::init(), libMesh::DiffSolver::initial_linear_tolerance, libMesh::INVALID_SOLVER_PACKAGE, libMesh::ExactSolution::l2_error(), libMesh::DiffSolver::max_linear_iterations, libMesh::DiffSolver::max_nonlinear_iterations, mesh, libMesh::out, libMesh::EquationSystems::print_info(), libMesh::MeshBase::print_info(), libMesh::DiffSolver::quiet, libMesh::DiffSolver::relative_residual_tolerance, libMesh::DiffSolver::relative_step_tolerance, libMesh::FEMSystem::solve(), libMesh::DifferentiableSystem::time_solver, libMesh::System::variable_number(), libMesh::DiffSolver::verbose, and libMesh::ExodusII_IO::write_equation_systems().

49 {
50  // Initialize libMesh.
51  LibMeshInit init (argc, argv);
52 
53  // This example requires a linear solver package.
54  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
55  "--enable-petsc, --enable-trilinos, or --enable-eigen");
56 
57  // Parse the input file
58  GetPot infile("vector_fe_ex4.in");
59 
60  // But allow the command line to override it.
61  infile.parse_command_line(argc, argv);
62 
63  // Read in parameters from the input file
64  const unsigned int grid_size = infile("grid_size", 2);
65 
66  // Skip higher-dimensional examples on a lower-dimensional libMesh build
67  libmesh_example_requires(3 <= LIBMESH_DIM, "3D support");
68 
69  // Create a mesh, with dimension to be overridden later, on the
70  // default MPI communicator.
71  Mesh mesh(init.comm());
72 
73  // Use the MeshTools::Generation mesh generator to create a uniform
74  // grid on the square [-1,1]^D. We must use at least TET10 or HEX20 elements
75  // for the Nedelec tetrahedral or hexahedral elements, respectively.
76  std::string elem_str =
77  command_line_value(std::string("element_type"),
78  std::string("HEX27"));
79 
80  // In general we expect O(h) convergence in *both* the L2 and
81  // H(curl) norms when using Nedelec elements. For more information
82  // on this topic, see the relevant results in other software [1-3]
83  // and the descriptions in [4,5]. In this particular example, we
84  // have observed O(h^2) convergence in L2 for HEX20/HEX27s, since the
85  // particular solution we seek is, just like the basis functions, constant
86  // along each of the cartesian axes. The basis functions are linear in the
87  // other two dimensions, improving the convergence speed.
88  //
89  // [1]: deal.ii, https://www.dealii.org/reports/nedelec/nedelec.pdf
90  // [2]: FEMPAR, https://www.sciencedirect.com/science/article/pii/S096599781831113X
91  // [3]: FEniCS, https://fenicsproject.org/pub/book/book/fenics-book-2011-06-24.pdf
92  // [4]: Monk, https://icerm.brown.edu/materials/Slides/tw-18-7/Finite_Element_Methods_for_Maxwells_Equations_%5D_Peter_Monk,_University_of_Delaware.pdf
93  // [5]: Hiptmair at al., https://www.sam.math.ethz.ch/sam_reports/reports_final/reports2009/2009-04_fp.pdf
94  libmesh_error_msg_if(elem_str != "TET10" && elem_str != "TET14" && elem_str != "HEX20" && elem_str != "HEX27",
95  "You entered: " << elem_str <<
96  " but this example must be run with TET10, TET14, HEX20 or HEX27.");
97 
99  grid_size,
100  grid_size,
101  grid_size,
102  -1., 1.,
103  -1., 1.,
104  -1., 1.,
105  Utility::string_to_enum<ElemType>(elem_str));
106 
107 
108  // Print information about the mesh to the screen.
109  mesh.print_info();
110 
111  // Create an equation systems object.
112  EquationSystems equation_systems (mesh);
113 
114  // Declare the system "CurlCurl" and its variables.
115  CurlCurlSystem & system =
116  equation_systems.add_system<CurlCurlSystem> ("CurlCurl");
117 
118  // This example only implements the steady-state problem
119  system.time_solver = std::make_unique<SteadySolver>(system);
120 
121  // Initialize the system
122  equation_systems.init();
123 
124  // And the nonlinear solver options
125  DiffSolver & solver = *(system.time_solver->diff_solver().get());
126  solver.quiet = infile("solver_quiet", true);
127  solver.verbose = !solver.quiet;
128  solver.max_nonlinear_iterations = infile("max_nonlinear_iterations", 15);
129  solver.relative_step_tolerance = infile("relative_step_tolerance", 1.e-3);
130  solver.relative_residual_tolerance = infile("relative_residual_tolerance", 1.0e-13);
131  solver.absolute_residual_tolerance = infile("absolute_residual_tolerance", 0.0);
132 
133  // And the linear solver options
134  solver.max_linear_iterations = infile("max_linear_iterations", 50000);
135  solver.initial_linear_tolerance = infile("initial_linear_tolerance", 1.e-10);
136 
137  // Print information about the system to the screen.
138  equation_systems.print_info();
139 
140  // Solve the system.
141  system.solve();
142 
143  ExactSolution exact_sol(equation_systems);
144 
145  SolutionFunction soln_func(system.variable_number("u"));
146  SolutionGradient soln_grad(system.variable_number("u"));
147 
148  // Build FunctionBase* containers to attach to the ExactSolution object.
149  std::vector<FunctionBase<Number> *> sols(1, &soln_func);
150  std::vector<FunctionBase<Gradient> *> grads(1, &soln_grad);
151 
152  exact_sol.attach_exact_values(sols);
153  exact_sol.attach_exact_derivs(grads);
154 
155  // Use higher quadrature order for more accurate error results
156  int extra_error_quadrature = infile("extra_error_quadrature", 2);
157  exact_sol.extra_quadrature_order(extra_error_quadrature);
158 
159  // Compute the error.
160  exact_sol.compute_error("CurlCurl", "u");
161 
162  // Print out the error values
163  libMesh::out << "L2-Error is: "
164  << exact_sol.l2_error("CurlCurl", "u")
165  << std::endl;
166  libMesh::out << "HCurl semi-norm error is: "
167  << exact_sol.error_norm("CurlCurl", "u", HCURL_SEMINORM)
168  << std::endl;
169  libMesh::out << "HCurl-Error is: "
170  << exact_sol.hcurl_error("CurlCurl", "u")
171  << std::endl;
172 
173 #ifdef LIBMESH_HAVE_EXODUS_API
174 
175  // We write the file in the ExodusII format.
176  ExodusII_IO(mesh).write_equation_systems("out.e", equation_systems);
177 
178 #endif // #ifdef LIBMESH_HAVE_EXODUS_API
179 
180  // All done.
181  return 0;
182 }
This class handles the computation of the L2 and/or H1 error for the Systems in the EquationSystems o...
This is the EquationSystems class.
bool quiet
The DiffSolver should not print anything to libMesh::out unless quiet is set to false; default is tru...
Definition: diff_solver.h:160
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
Definition: exodusII_io.h:52
unsigned int max_nonlinear_iterations
The DiffSolver should exit in failure if max_nonlinear_iterations is exceeded and continue_after_max_...
Definition: diff_solver.h:154
Real absolute_residual_tolerance
The DiffSolver should exit after the residual is reduced to either less than absolute_residual_tolera...
Definition: diff_solver.h:189
std::unique_ptr< TimeSolver > time_solver
A pointer to the solver object we&#39;re going to use.
Definition: diff_system.h:253
MeshBase & mesh
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
Definition: libmesh.h:90
unsigned int variable_number(std::string_view var) const
Definition: system.C:1609
SolverPackage default_solver_package()
Definition: libmesh.C:1117
This is a generic class that defines a solver to handle ImplicitSystem classes, including NonlinearIm...
Definition: diff_solver.h:68
T command_line_value(const std::string &, T)
Definition: libmesh.C:1024
unsigned int max_linear_iterations
Each linear solver step should exit after max_linear_iterations is exceeded.
Definition: diff_solver.h:146
virtual void write_equation_systems(const std::string &fname, const EquationSystems &es, const std::set< std::string > *system_names=nullptr) override
Writes out the solution for no specific time or timestep.
Definition: exodusII_io.C:2033
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
void init(triangulateio &t)
Initializes the fields of t to nullptr/0 as necessary.
FEMSystem, TimeSolver and NewtonSolver will handle most tasks, but we must specify element residuals...
bool verbose
The DiffSolver may print a lot more to libMesh::out if verbose is set to true; default is false...
Definition: diff_solver.h:166
OStreamProxy out
virtual void solve() override
Invokes the solver associated with the system.
Definition: fem_system.C:1076
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
double initial_linear_tolerance
Any required linear solves will at first be done with this tolerance; the DiffSolver may tighten the ...
Definition: diff_solver.h:208
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.
Real relative_residual_tolerance
Definition: diff_solver.h:190