35 #include "libmesh/libmesh.h" 36 #include "libmesh/mesh.h" 37 #include "libmesh/mesh_generation.h" 38 #include "libmesh/exodusII_io.h" 39 #include "libmesh/gnuplot_io.h" 40 #include "libmesh/linear_implicit_system.h" 41 #include "libmesh/equation_systems.h" 44 #include "libmesh/fe.h" 47 #include "libmesh/quadrature_gauss.h" 51 #include "libmesh/dof_map.h" 55 #include "libmesh/sparse_matrix.h" 56 #include "libmesh/numeric_vector.h" 57 #include "libmesh/dense_matrix.h" 58 #include "libmesh/dense_vector.h" 63 #include "libmesh/perf_log.h" 66 #include "libmesh/elem.h" 68 #include "libmesh/string_to_enum.h" 69 #include "libmesh/getpot.h" 70 #include "libmesh/enum_solver_package.h" 84 const std::string & system_name);
92 int main (
int argc,
char ** argv)
99 "--enable-petsc, --enable-trilinos, or --enable-eigen");
105 GetPot command_line (argc, argv);
108 libmesh_error_msg_if(argc < 3,
"Usage:\n" <<
"\t " << argv[0] <<
" -d 2(3)" <<
" -n 15");
114 for (
int i=1; i<argc; i++)
125 libmesh_example_requires(
dim <= LIBMESH_DIM,
"2D/3D support");
135 std::string order =
"SECOND";
140 std::string family =
"LAGRANGE";
145 libmesh_error_msg_if(((family ==
"MONOMIAL") || (family ==
"XYZ")) &&
147 "This example requires a C^0 (or higher) FE basis.");
155 Real halfwidth =
dim > 1 ? 1. : 0.;
156 Real halfheight =
dim > 2 ? 1. : 0.;
158 if ((family ==
"LAGRANGE") && (order ==
"FIRST"))
167 -halfwidth, halfwidth,
168 -halfheight, halfheight,
180 -halfwidth, halfwidth,
181 -halfheight, halfheight,
186 for (
auto & elem :
mesh.element_ptr_range())
188 const Point cent = elem->vertex_average();
191 if ((cent(0) > 0) == (cent(1) > 0))
192 elem->subdomain_id() = 1;
194 else if (cent(0) > 0)
195 elem->subdomain_id() = 1;
210 std::set<subdomain_id_type> active_subdomains;
215 active_subdomains.
clear(); active_subdomains.insert(0);
217 Utility::string_to_enum<Order> (order),
218 Utility::string_to_enum<FEFamily>(family),
223 active_subdomains.clear(); active_subdomains.insert(1);
225 Utility::string_to_enum<Order> (order),
226 Utility::string_to_enum<FEFamily>(family),
234 equation_systems.
init();
241 equation_systems.
get_system(
"Poisson").solve();
252 #ifdef LIBMESH_HAVE_EXODUS_API 254 "out_3.e" :
"out_2.e", equation_systems);
255 #endif // #ifdef LIBMESH_HAVE_EXODUS_API 274 const std::string & libmesh_dbg_var(system_name))
278 libmesh_assert_equal_to (system_name,
"Poisson");
284 PerfLog perf_log (
"Matrix Assembly");
303 FEType fe_type = dof_map.variable_type(0);
315 fe->attach_quadrature_rule (&qrule);
328 fe_face->attach_quadrature_rule (&qface);
334 const std::vector<Real> & JxW = fe->get_JxW();
339 const std::vector<Point> & q_point = fe->get_xyz();
342 const std::vector<std::vector<Real>> & phi = fe->get_phi();
346 const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
358 std::vector<dof_id_type> dof_indices, dof_indices2;
371 for (
const auto & elem :
as_range(
mesh.local_elements_begin(),
372 mesh.local_elements_end()))
377 perf_log.
push(
"elem init");
383 dof_map.dof_indices (elem, dof_indices, 0);
384 dof_map.dof_indices (elem, dof_indices2, 1);
404 Ke.
resize (std::max(dof_indices.size(), dof_indices2.size()),
405 std::max(dof_indices.size(), dof_indices2.size()));
407 Fe.
resize (std::max(dof_indices.size(), dof_indices2.size()));
412 perf_log.
pop(
"elem init");
423 perf_log.
push (
"Ke");
425 for (
unsigned int qp=0; qp<qrule.
n_points(); qp++)
426 for (std::size_t i=0; i<phi.size(); i++)
427 for (std::size_t j=0; j<phi.size(); j++)
428 Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
439 perf_log.
push (
"Fe");
441 for (
unsigned int qp=0; qp<qrule.
n_points(); qp++)
457 const Real x = q_point[qp](0);
459 const Real y = q_point[qp](1);
464 const Real z = q_point[qp](2);
468 const Real eps = 1.e-3;
488 fxy = (0.25*
pi*
pi)*sin(.5*
pi*x);
492 fxy = - (uxx + uyy + ((
dim==2) ? 0. : uzz));
496 for (std::size_t i=0; i<phi.size(); i++)
497 Fe(i) += JxW[qp]*fxy*phi[i][qp];
512 LOG_SCOPE_WITH(
"BCs",
"", perf_log);
517 for (
auto side : elem->side_index_range())
518 if ((elem->neighbor_ptr(side) ==
nullptr) ||
519 (elem->neighbor_ptr(side)->subdomain_id() != elem->subdomain_id()))
524 const Real penalty = 1.e10;
528 const std::vector<std::vector<Real>> & phi_face = fe_face->get_phi();
532 const std::vector<Real> & JxW_face = fe_face->get_JxW();
537 const std::vector<Point> & qface_point = fe_face->get_xyz();
541 fe_face->reinit(elem, side);
544 for (
unsigned int qp=0; qp<qface.
n_points(); qp++)
548 const Real xf = qface_point[qp](0);
550 const Real yf = qface_point[qp](1);
555 const Real zf = qface_point[qp](2);
565 for (std::size_t i=0; i<phi_face.size(); i++)
566 for (std::size_t j=0; j<phi_face.size(); j++)
567 Ke(i,j) += JxW_face[qp]*penalty*phi_face[i][qp]*phi_face[j][qp];
571 for (std::size_t i=0; i<phi_face.size(); i++)
572 Fe(i) += JxW_face[qp]*penalty*
value*phi_face[i][qp];
584 LOG_SCOPE_WITH(
"matrix insertion",
"", perf_log);
586 if (dof_indices.size())
592 if (dof_indices2.size())
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
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.
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.
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
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.
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 print_info(std::ostream &os=libMesh::out, const unsigned int verbosity=0, const bool global=true) const
Prints relevant information about the mesh.
SimpleRange< IndexType > as_range(const std::pair< IndexType, IndexType > &p)
Helper function that allows us to treat a homogenous pair as a range.
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
int main(int argc, char **argv)
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
virtual void clear() override
Clear all the data structures associated with 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.
Real exact_solution(const Real x, const Real y=0., const Real z=0.)
This is the exact solution that we are trying to obtain.
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
processor_id_type processor_id() const
const DofMap & get_dof_map() const
A Point defines a location in LIBMESH_DIM dimensional Real space.
const SparseMatrix< Number > & get_system_matrix() const
void assemble_poisson(EquationSystems &es, const std::string &system_name)