46 #include "libmesh/libmesh.h" 47 #include "libmesh/mesh.h" 48 #include "libmesh/mesh_generation.h" 49 #include "libmesh/exodusII_io.h" 50 #include "libmesh/gnuplot_io.h" 51 #include "libmesh/linear_implicit_system.h" 52 #include "libmesh/equation_systems.h" 55 #include "libmesh/fe.h" 58 #include "libmesh/quadrature_gauss.h" 62 #include "libmesh/dof_map.h" 66 #include "libmesh/sparse_matrix.h" 67 #include "libmesh/numeric_vector.h" 68 #include "libmesh/dense_matrix.h" 69 #include "libmesh/dense_vector.h" 74 #include "libmesh/perf_log.h" 77 #include "libmesh/elem.h" 80 #include "libmesh/dirichlet_boundaries.h" 81 #include "libmesh/analytic_function.h" 83 #include "libmesh/string_to_enum.h" 84 #include "libmesh/getpot.h" 85 #include "libmesh/enum_solver_package.h" 99 const std::string & system_name);
112 (LIBMESH_DIM>1)?p(1):0,
113 (LIBMESH_DIM>2)?p(2):0);
117 int main (
int argc,
char ** argv)
124 "--enable-petsc, --enable-trilinos, or --enable-eigen");
130 libmesh_error_msg_if(argc < 3,
"Usage:\n" <<
"\t " << argv[0] <<
" -d 2(3)" <<
" -n 15");
136 for (
int i=1; i<argc; i++)
147 libmesh_example_requires(
dim <= LIBMESH_DIM,
"2D/3D support");
150 #ifndef LIBMESH_ENABLE_DIRICHLET 151 libmesh_example_requires(
false,
"--enable-dirichlet");
159 std::string order =
"SECOND";
164 std::string family =
"LAGRANGE";
169 libmesh_error_msg_if((family ==
"MONOMIAL") || (family ==
"XYZ"),
170 "ex4 currently requires a C^0 (or higher) FE basis.");
182 Real halfwidth =
dim > 1 ? 1. : 0.;
183 Real halfheight =
dim > 2 ? 1. : 0.;
185 if ((family ==
"LAGRANGE") && (order ==
"FIRST"))
194 -halfwidth, halfwidth,
195 -halfheight, halfheight,
207 -halfwidth, halfwidth,
208 -halfheight, halfheight,
230 Utility::string_to_enum<Order> (order),
231 Utility::string_to_enum<FEFamily>(family));
242 std::set<boundary_id_type> boundary_ids;
244 boundary_ids.insert(0);
245 boundary_ids.insert(1);
249 boundary_ids.insert(2);
250 boundary_ids.insert(3);
255 boundary_ids.insert(4);
256 boundary_ids.insert(5);
263 #ifdef LIBMESH_ENABLE_DIRICHLET 270 (boundary_ids, {u_var}, exact_solution_object);
278 equation_systems.
init();
294 #ifdef LIBMESH_HAVE_EXODUS_API 298 "out_3.e" :
"out_2.e", equation_systems);
300 #endif // #ifdef LIBMESH_HAVE_EXODUS_API 314 const std::string & libmesh_dbg_var(system_name))
318 libmesh_assert_equal_to (system_name,
"Poisson");
324 PerfLog perf_log (
"Matrix Assembly");
343 FEType fe_type = dof_map.variable_type(0);
355 fe->attach_quadrature_rule (&qrule);
368 fe_face->attach_quadrature_rule (&qface);
374 const std::vector<Real> & JxW = fe->get_JxW();
379 const std::vector<Point> & q_point = fe->get_xyz();
382 const std::vector<std::vector<Real>> & phi = fe->get_phi();
386 const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
398 std::vector<dof_id_type> dof_indices;
407 for (
const auto & elem :
mesh.active_local_element_ptr_range())
412 perf_log.
push(
"elem init");
418 dof_map.dof_indices (elem, dof_indices);
425 const unsigned int n_dofs =
426 cast_int<unsigned int>(dof_indices.size());
436 libmesh_assert_equal_to (n_dofs, phi.size());
444 Ke.
resize (n_dofs, n_dofs);
451 perf_log.
pop(
"elem init");
462 perf_log.
push (
"Ke");
464 for (
unsigned int qp=0; qp<qrule.
n_points(); qp++)
465 for (
unsigned int i=0; i != n_dofs; i++)
466 for (
unsigned int j=0; j != n_dofs; j++)
467 Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
478 perf_log.
push (
"Fe");
480 for (
unsigned int qp=0; qp<qrule.
n_points(); qp++)
496 const Real x = q_point[qp](0);
498 const Real y = q_point[qp](1);
503 const Real z = q_point[qp](2);
507 const Real eps = 1.e-3;
527 fxy = (0.25*
pi*
pi)*sin(.5*
pi*x);
531 fxy = - (uxx + uyy + ((
dim==2) ? 0. : uzz));
535 for (
unsigned int i=0; i != n_dofs; i++)
536 Fe(i) += JxW[qp]*fxy*phi[i][qp];
546 dof_map.heterogenously_constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
554 LOG_SCOPE_WITH(
"matrix insertion",
"", perf_log);
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
int main(int argc, char **argv)
T command_line_next(std::string name, T default_value)
Use GetPot's search()/next() functions to get following arguments from the command line...
This is the EquationSystems class.
void pop(const char *label, const char *header="")
Pop the event label off the stack, resuming any lower event.
void exact_solution_wrapper(DenseVector< Number > &output, const Point &p, const Real)
virtual void write_equation_systems(const std::string &, const EquationSystems &, const std::set< std::string > *system_names=nullptr)
This method implements writing a mesh with data to a specified file where the data is taken from the ...
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
Manages consistently variables, degrees of freedom, coefficient vectors, matrices and linear solvers ...
void resize(const unsigned int n)
Resize the vector.
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.
NumericVector< Number > * rhs
The system matrix.
This class allows one to associate Dirichlet boundary values with a given set of mesh boundary ids an...
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
The libMesh namespace provides an interface to certain functionality in the library.
virtual void solve() override
Assembles & solves the linear system A*x=b.
const T_sys & get_system(std::string_view name) const
Wraps a function pointer into a FunctionBase object.
This is the MeshBase class.
SolverPackage default_solver_package()
The PerfLog class allows monitoring of specific events.
This class handles the numbering of degrees of freedom on a mesh.
This class implements writing meshes using GNUplot, designed for use only with 1D meshes...
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.
void assemble_poisson(EquationSystems &es, const std::string &system_name)
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.
void print_info(std::ostream &os=libMesh::out, const unsigned int verbosity=0, const bool global=true) const
Prints relevant information about the mesh.
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.
unsigned int add_variable(std::string_view var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=nullptr)
Adds the variable var to the list of variables for this system.
void push(const char *label, const char *header="")
Push the event label onto the stack, pausing any active event.
unsigned int n_points() const
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.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
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().
This class implements specific orders of Gauss quadrature.
unsigned int mesh_dimension() const
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 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.
const DofMap & get_dof_map() const
A Point defines a location in LIBMESH_DIM dimensional Real space.
const SparseMatrix< Number > & get_system_matrix() const
Real exact_solution(const Real x, const Real y, const Real z=0.)
This is the exact solution that we are trying to obtain.