libMesh
Functions
miscellaneous_ex1.C File Reference

Go to the source code of this file.

Functions

void assemble_wave (EquationSystems &es, const std::string &system_name)
 
int main (int argc, char **argv)
 
void assemble_wave (EquationSystems &es, const std::string &libmesh_dbg_var(system_name))
 

Function Documentation

void assemble_wave ( EquationSystems es,
const std::string &  system_name 
)

Definition at line 322 of file transient_ex2.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::DofMap::dof_indices(), libMesh::Parameters::get(), libMesh::EquationSystems::get_mesh(), libMesh::EquationSystems::get_system(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::QBase::n_points(), libMesh::EquationSystems::parameters, libMesh::Real, libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), libMesh::SECOND, and side.

Referenced by main().

324 {
325  // It is a good idea to make sure we are assembling
326  // the proper system.
327  libmesh_assert_equal_to (system_name, "Wave");
328 
329  // Get a constant reference to the mesh object.
330  const MeshBase & mesh = es.get_mesh();
331 
332  // The dimension that we are running.
333  const unsigned int dim = mesh.mesh_dimension();
334 
335  // Copy the speed of sound and fluid density
336  // to a local variable.
337  const Real speed = es.parameters.get<Real>("speed");
338  const Real rho = es.parameters.get<Real>("fluid density");
339 
340  // Get a reference to our system, as before.
341  NewmarkSystem & t_system = es.get_system<NewmarkSystem> (system_name);
342 
343  // Get a constant reference to the Finite Element type
344  // for the first (and only) variable in the system.
345  FEType fe_type = t_system.get_dof_map().variable_type(0);
346 
347  // In here, we will add the element matrices to the
348  // @e additional matrices "stiffness_mass" and "damping"
349  // and the additional vector "force", not to the members
350  // "matrix" and "rhs". Therefore, get writable
351  // references to them.
352  SparseMatrix<Number> & stiffness = t_system.get_matrix("stiffness");
353  SparseMatrix<Number> & damping = t_system.get_matrix("damping");
354  SparseMatrix<Number> & mass = t_system.get_matrix("mass");
355  NumericVector<Number> & force = t_system.get_vector("force");
356 
357  // Some solver packages (PETSc) are especially picky about
358  // allocating sparsity structure and truly assigning values
359  // to this structure. Namely, matrix additions, as performed
360  // later, exhibit acceptable performance only for identical
361  // sparsity structures. Therefore, explicitly zero the
362  // values in the collective matrix, so that matrix additions
363  // encounter identical sparsity structures.
364  SparseMatrix<Number> & matrix = *t_system.matrix;
365  DenseMatrix<Number> zero_matrix;
366 
367  // Build a Finite Element object of the specified type. Since the
368  // FEBase::build() member dynamically creates memory we will
369  // store the object as a UniquePtr<FEBase>. This can be thought
370  // of as a pointer that will clean up after itself.
371  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
372 
373  // A 2nd order Gauss quadrature rule for numerical integration.
374  QGauss qrule (dim, SECOND);
375 
376  // Tell the finite element object to use our quadrature rule.
377  fe->attach_quadrature_rule (&qrule);
378 
379  // The element Jacobian * quadrature weight at each integration point.
380  const std::vector<Real> & JxW = fe->get_JxW();
381 
382  // The element shape functions evaluated at the quadrature points.
383  const std::vector<std::vector<Real>> & phi = fe->get_phi();
384 
385  // The element shape function gradients evaluated at the quadrature
386  // points.
387  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
388 
389  // A reference to the DofMap object for this system. The DofMap
390  // object handles the index translation from node and element numbers
391  // to degree of freedom numbers.
392  const DofMap & dof_map = t_system.get_dof_map();
393 
394  // The element mass, damping and stiffness matrices
395  // and the element contribution to the rhs.
396  DenseMatrix<Number> Ke, Ce, Me;
398 
399  // This vector will hold the degree of freedom indices for
400  // the element. These define where in the global system
401  // the element degrees of freedom get mapped.
402  std::vector<dof_id_type> dof_indices;
403 
404  // Now we will loop over all the elements in the mesh.
405  // We will compute the element matrix and right-hand-side
406  // contribution.
407  for (const auto & elem : mesh.active_local_element_ptr_range())
408  {
409  // Get the degree of freedom indices for the
410  // current element. These define where in the global
411  // matrix and right-hand-side this element will
412  // contribute to.
413  dof_map.dof_indices (elem, dof_indices);
414 
415  // Compute the element-specific data for the current
416  // element. This involves computing the location of the
417  // quadrature points (q_point) and the shape functions
418  // (phi, dphi) for the current element.
419  fe->reinit (elem);
420 
421  // Zero the element matrices and rhs before
422  // summing them. We use the resize member here because
423  // the number of degrees of freedom might have changed from
424  // the last element. Note that this will be the case if the
425  // element type is different (i.e. the last element was HEX8
426  // and now have a PRISM6).
427  {
428  const unsigned int n_dof_indices = dof_indices.size();
429 
430  Ke.resize (n_dof_indices, n_dof_indices);
431  Ce.resize (n_dof_indices, n_dof_indices);
432  Me.resize (n_dof_indices, n_dof_indices);
433  zero_matrix.resize (n_dof_indices, n_dof_indices);
434  Fe.resize (n_dof_indices);
435  }
436 
437  // Now loop over the quadrature points. This handles
438  // the numeric integration.
439  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
440  {
441  // Now we will build the element matrix. This involves
442  // a double loop to integrate the test functions (i) against
443  // the trial functions (j).
444  for (std::size_t i=0; i<phi.size(); i++)
445  for (std::size_t j=0; j<phi.size(); j++)
446  {
447  Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
448  Me(i,j) += JxW[qp]*phi[i][qp]*phi[j][qp]
449  *1./(speed*speed);
450  } // end of the matrix summation loop
451  } // end of quadrature point loop
452 
453  // Now compute the contribution to the element matrix and the
454  // right-hand-side vector if the current element lies on the
455  // boundary.
456  {
457  // In this example no natural boundary conditions will
458  // be considered. The code is left here so it can easily
459  // be extended.
460  //
461  // don't do this for any side
462  for (auto side : elem->side_index_range())
463  if (!true)
464  // if (elem->neighbor_ptr(side) == libmesh_nullptr)
465  {
466  // Declare a special finite element object for
467  // boundary integration.
468  UniquePtr<FEBase> fe_face (FEBase::build(dim, fe_type));
469 
470  // Boundary integration requires one quadrature rule,
471  // with dimensionality one less than the dimensionality
472  // of the element.
473  QGauss qface(dim-1, SECOND);
474 
475  // Tell the finite element object to use our
476  // quadrature rule.
477  fe_face->attach_quadrature_rule (&qface);
478 
479  // The value of the shape functions at the quadrature
480  // points.
481  const std::vector<std::vector<Real>> & phi_face = fe_face->get_phi();
482 
483  // The Jacobian * Quadrature Weight at the quadrature
484  // points on the face.
485  const std::vector<Real> & JxW_face = fe_face->get_JxW();
486 
487  // Compute the shape function values on the element
488  // face.
489  fe_face->reinit(elem, side);
490 
491  // Here we consider a normal acceleration acc_n=1 applied to
492  // the whole boundary of our mesh.
493  const Real acc_n_value = 1.0;
494 
495  // Loop over the face quadrature points for integration.
496  for (unsigned int qp=0; qp<qface.n_points(); qp++)
497  {
498  // Right-hand-side contribution due to prescribed
499  // normal acceleration.
500  for (std::size_t i=0; i<phi_face.size(); i++)
501  {
502  Fe(i) += acc_n_value*rho
503  *phi_face[i][qp]*JxW_face[qp];
504  }
505  } // end face quadrature point loop
506  } // end if (elem->neighbor_ptr(side) == libmesh_nullptr)
507 
508  // In this example the Dirichlet boundary conditions will be
509  // imposed via penalty method after the
510  // system is assembled.
511 
512  } // end boundary condition section
513 
514  // If this assembly program were to be used on an adaptive mesh,
515  // we would have to apply any hanging node constraint equations
516  // by uncommenting the following lines:
517  // std::vector<unsigned int> dof_indicesC = dof_indices;
518  // std::vector<unsigned int> dof_indicesM = dof_indices;
519  // dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
520  // dof_map.constrain_element_matrix (Ce, dof_indicesC);
521  // dof_map.constrain_element_matrix (Me, dof_indicesM);
522 
523  // Finally, simply add the contributions to the additional
524  // matrices and vector.
525  stiffness.add_matrix (Ke, dof_indices);
526  damping.add_matrix (Ce, dof_indices);
527  mass.add_matrix (Me, dof_indices);
528 
529  force.add_vector (Fe, dof_indices);
530 
531  // For the overall matrix, explicitly zero the entries where
532  // we added values in the other ones, so that we have
533  // identical sparsity footprints.
534  matrix.add_matrix(zero_matrix, dof_indices);
535 
536  } // end of element loop
537 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
unsigned int dim
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]...
unsigned short int side
Definition: xdr_io.C:49
MeshBase & mesh
const T & get(const std::string &) const
Definition: parameters.h:430
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.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
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 contains a specific system class.
This class implements specific orders of Gauss quadrature.
Parameters parameters
Data structure holding arbitrary parameters.
const MeshBase & get_mesh() const
const T_sys & get_system(const std::string &name) const
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Fills the vector di with the global degree of freedom indices for the element.
Definition: dof_map.C:1917
void assemble_wave ( EquationSystems es,
const std::string &  libmesh_dbg_varsystem_name 
)

Definition at line 222 of file miscellaneous_ex1.C.

References std::abs(), libMesh::MeshBase::active_local_element_ptr_range(), libMesh::DenseMatrix< T >::add(), libMesh::NumericVector< T >::add(), libMesh::SparseMatrix< T >::add_matrix(), libMesh::NumericVector< T >::add_vector(), libMesh::FEGenericBase< OutputType >::build(), libMesh::FEGenericBase< OutputType >::build_InfFE(), dim, libMesh::Parameters::get(), libMesh::System::get_dof_map(), libMesh::FEGenericBase< OutputType >::get_dphase(), libMesh::FEGenericBase< OutputType >::get_dphi(), libMesh::FEAbstract::get_JxW(), libMesh::EquationSystems::get_mesh(), libMesh::FEGenericBase< OutputType >::get_phi(), libMesh::FEGenericBase< OutputType >::get_Sobolev_dweight(), libMesh::FEGenericBase< OutputType >::get_Sobolev_weight(), libMesh::EquationSystems::get_system(), libMesh::libmesh_ignore(), libmesh_nullptr, libMesh::MeshBase::local_node_ptr_range(), libMesh::ImplicitSystem::matrix, mesh, libMesh::MeshBase::mesh_dimension(), libMesh::FEAbstract::n_quadrature_points(), libMesh::FEAbstract::n_shape_functions(), libMesh::EquationSystems::parameters, libMesh::Real, libMesh::FEAbstract::reinit(), libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), libMesh::ExplicitSystem::rhs, libMesh::SECOND, libMesh::TOLERANCE, and libMesh::MeshTools::weight().

224 {
225  // It is a good idea to make sure we are assembling
226  // the proper system.
227  libmesh_assert_equal_to (system_name, "Wave");
228 
229  // Avoid unused variable warnings when compiling without infinite
230  // elements enabled.
231  libmesh_ignore(es);
232 
233 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
234 
235  // Get a constant reference to the mesh object.
236  const MeshBase & mesh = es.get_mesh();
237 
238  // Get a reference to the system we are solving.
240 
241  // A reference to the DofMap object for this system. The DofMap
242  // object handles the index translation from node and element numbers
243  // to degree of freedom numbers.
244  const DofMap & dof_map = system.get_dof_map();
245 
246  // The dimension that we are running.
247  const unsigned int dim = mesh.mesh_dimension();
248 
249  // Copy the speed of sound to a local variable.
250  const Real speed = es.parameters.get<Real>("speed");
251 
252  // Get a constant reference to the Finite Element type
253  // for the first (and only) variable in the system.
254  const FEType & fe_type = dof_map.variable_type(0);
255 
256  // Build a Finite Element object of the specified type. Since the
257  // FEBase::build() member dynamically creates memory we will
258  // store the object as a UniquePtr<FEBase>.
259  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
260 
261  // Do the same for an infinite element.
262  UniquePtr<FEBase> inf_fe (FEBase::build_InfFE(dim, fe_type));
263 
264  // A 2nd order Gauss quadrature rule for numerical integration.
265  QGauss qrule (dim, SECOND);
266 
267  // Tell the finite element object to use our quadrature rule.
268  fe->attach_quadrature_rule (&qrule);
269 
270  // Due to its internal structure, the infinite element handles
271  // quadrature rules differently. It takes the quadrature
272  // rule which has been initialized for the FE object, but
273  // creates suitable quadrature rules by @e itself. The user
274  // need not worry about this.
275  inf_fe->attach_quadrature_rule (&qrule);
276 
277  // Define data structures to contain the element matrix
278  // and right-hand-side vector contribution. Following
279  // basic finite element terminology we will denote these
280  // "Ke", "Ce", "Me", and "Fe" for the stiffness, damping
281  // and mass matrices, and the load vector. Note that in
282  // Acoustics, these descriptors though do @e not match the
283  // true physical meaning of the projectors. The final
284  // overall system, however, resembles the conventional
285  // notation again.
290 
291  // This vector will hold the degree of freedom indices for
292  // the element. These define where in the global system
293  // the element degrees of freedom get mapped.
294  std::vector<dof_id_type> dof_indices;
295 
296  // Now we will loop over all the elements in the mesh.
297  // We will compute the element matrix and right-hand-side
298  // contribution.
299  for (const auto & elem : mesh.active_local_element_ptr_range())
300  {
301  // Get the degree of freedom indices for the
302  // current element. These define where in the global
303  // matrix and right-hand-side this element will
304  // contribute to.
305  dof_map.dof_indices (elem, dof_indices);
306 
307  // The mesh contains both finite and infinite elements. These
308  // elements are handled through different classes, namely
309  // FE and InfFE, respectively. However, since both
310  // are derived from FEBase, they share the same interface,
311  // and overall burden of coding is @e greatly reduced through
312  // using a pointer, which is adjusted appropriately to the
313  // current element type.
314  FEBase * cfe = libmesh_nullptr;
315 
316  // This here is almost the only place where we need to
317  // distinguish between finite and infinite elements.
318  // For faster computation, however, different approaches
319  // may be feasible.
320  //
321  // Up to now, we do not know what kind of element we
322  // have. Aske the element of what type it is:
323  if (elem->infinite())
324  {
325  // We have an infinite element. Let cfe point
326  // to our InfFE object. This is handled through
327  // a UniquePtr. Through the UniquePtr::get() we "borrow"
328  // the pointer, while the UniquePtr inf_fe is
329  // still in charge of memory management.
330  cfe = inf_fe.get();
331  }
332  else
333  {
334  // This is a conventional finite element. Let fe handle it.
335  cfe = fe.get();
336 
337  // Boundary conditions.
338  // Here we just zero the rhs-vector. For natural boundary
339  // conditions check e.g. previous examples.
340  {
341  // Zero the RHS for this element.
342  Fe.resize (dof_indices.size());
343 
344  system.rhs->add_vector (Fe, dof_indices);
345  } // end boundary condition section
346  } // else if (elem->infinite())
347 
348  // This is slightly different from the Poisson solver:
349  // Since the finite element object may change, we have to
350  // initialize the constant references to the data fields
351  // each time again, when a new element is processed.
352  //
353  // The element Jacobian * quadrature weight at each integration point.
354  const std::vector<Real> & JxW = cfe->get_JxW();
355 
356  // The element shape functions evaluated at the quadrature points.
357  const std::vector<std::vector<Real>> & phi = cfe->get_phi();
358 
359  // The element shape function gradients evaluated at the quadrature
360  // points.
361  const std::vector<std::vector<RealGradient>> & dphi = cfe->get_dphi();
362 
363  // The infinite elements need more data fields than conventional FE.
364  // These are the gradients of the phase term dphase, an additional
365  // radial weight for the test functions Sobolev_weight, and its
366  // gradient.
367  //
368  // Note that these data fields are also initialized appropriately by
369  // the FE method, so that the weak form (below) is valid for @e both
370  // finite and infinite elements.
371  const std::vector<RealGradient> & dphase = cfe->get_dphase();
372  const std::vector<Real> & weight = cfe->get_Sobolev_weight();
373  const std::vector<RealGradient> & dweight = cfe->get_Sobolev_dweight();
374 
375  // Now this is all independent of whether we use an FE
376  // or an InfFE. Nice, hm? ;-)
377  //
378  // Compute the element-specific data, as described
379  // in previous examples.
380  cfe->reinit (elem);
381 
382  // Zero the element matrices. Boundary conditions were already
383  // processed in the FE-only section, see above.
384  Ke.resize (dof_indices.size(), dof_indices.size());
385  Ce.resize (dof_indices.size(), dof_indices.size());
386  Me.resize (dof_indices.size(), dof_indices.size());
387 
388  // The total number of quadrature points for infinite elements
389  // @e has to be determined in a different way, compared to
390  // conventional finite elements. This type of access is also
391  // valid for finite elements, so this can safely be used
392  // anytime, instead of asking the quadrature rule, as
393  // seen in previous examples.
394  unsigned int max_qp = cfe->n_quadrature_points();
395 
396  // Loop over the quadrature points.
397  for (unsigned int qp=0; qp<max_qp; qp++)
398  {
399  // Similar to the modified access to the number of quadrature
400  // points, the number of shape functions may also be obtained
401  // in a different manner. This offers the great advantage
402  // of being valid for both finite and infinite elements.
403  const unsigned int n_sf = cfe->n_shape_functions();
404 
405  // Now we will build the element matrices. Since the infinite
406  // elements are based on a Petrov-Galerkin scheme, the
407  // resulting system matrices are non-symmetric. The additional
408  // weight, described before, is part of the trial space.
409  //
410  // For the finite elements, though, these matrices are symmetric
411  // just as we know them, since the additional fields dphase,
412  // weight, and dweight are initialized appropriately.
413  //
414  // test functions: weight[qp]*phi[i][qp]
415  // trial functions: phi[j][qp]
416  // phase term: phase[qp]
417  //
418  // derivatives are similar, but note that these are of type
419  // Point, not of type Real.
420  for (unsigned int i=0; i<n_sf; i++)
421  for (unsigned int j=0; j<n_sf; j++)
422  {
423  // (ndt*Ht + nHt*d) * nH
424  Ke(i,j) +=
425  (
426  (dweight[qp] * phi[i][qp] // Point * Real = Point
427  + // +
428  dphi[i][qp] * weight[qp] // Point * Real = Point
429  ) * dphi[j][qp]
430  ) * JxW[qp];
431 
432  // (d*Ht*nmut*nH - ndt*nmu*Ht*H - d*nHt*nmu*H)
433  Ce(i,j) +=
434  (
435  (dphase[qp] * dphi[j][qp]) // (Point * Point) = Real
436  * weight[qp] * phi[i][qp] // * Real * Real = Real
437  - // -
438  (dweight[qp] * dphase[qp]) // (Point * Point) = Real
439  * phi[i][qp] * phi[j][qp] // * Real * Real = Real
440  - // -
441  (dphi[i][qp] * dphase[qp]) // (Point * Point) = Real
442  * weight[qp] * phi[j][qp] // * Real * Real = Real
443  ) * JxW[qp];
444 
445  // (d*Ht*H * (1 - nmut*nmu))
446  Me(i,j) +=
447  (
448  (1. - (dphase[qp] * dphase[qp])) // (Real - (Point * Point)) = Real
449  * phi[i][qp] * phi[j][qp] * weight[qp] // * Real * Real * Real = Real
450  ) * JxW[qp];
451 
452  } // end of the matrix summation loop
453  } // end of quadrature point loop
454 
455  // The element matrices are now built for this element.
456  // Collect them in Ke, and then add them to the global matrix.
457  // The SparseMatrix::add_matrix() member does this for us.
458  Ke.add(1./speed , Ce);
459  Ke.add(1./(speed*speed), Me);
460 
461  // If this assembly program were to be used on an adaptive mesh,
462  // we would have to apply any hanging node constraint equations
463  dof_map.constrain_element_matrix(Ke, dof_indices);
464 
465  system.matrix->add_matrix (Ke, dof_indices);
466  } // end of element loop
467 
468  // Note that we have not applied any boundary conditions so far.
469  // Here we apply a unit load at the node located at (0,0,0).
470  for (const auto & node : mesh.local_node_ptr_range())
471  if (std::abs((*node)(0)) < TOLERANCE &&
472  std::abs((*node)(1)) < TOLERANCE &&
473  std::abs((*node)(2)) < TOLERANCE)
474  {
475  // The global number of the respective degree of freedom.
476  unsigned int dn = node->dof_number(0,0,0);
477 
478  system.rhs->add (dn, 1.);
479  }
480 
481 #else
482 
483  // dummy assert
484  libmesh_assert_not_equal_to (es.get_mesh().mesh_dimension(), 1);
485 
486 #endif //ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
487 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
double abs(double a)
virtual SimpleRange< node_iterator > local_node_ptr_range()=0
unsigned int dim
virtual unsigned int n_quadrature_points() const =0
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
const class libmesh_nullptr_t libmesh_nullptr
NumericVector< Number > * rhs
The system matrix.
const T & get(const std::string &) const
Definition: parameters.h:430
static const Real TOLERANCE
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 unsigned int n_shape_functions() const =0
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
const std::vector< Real > & get_JxW() const
Definition: fe_abstract.h:237
const std::vector< Real > & get_Sobolev_weight() const
Definition: fe_base.h:420
dof_id_type weight(const MeshBase &mesh, const processor_id_type pid)
Definition: mesh_tools.C:245
boostcopy::enable_if_c< ScalarTraits< T2 >::value, void >::type add(const T2 factor, const DenseMatrix< T3 > &mat)
Adds factor times mat to this matrix.
Definition: dense_matrix.h:883
const std::vector< std::vector< OutputShape > > & get_phi() const
Definition: fe_base.h:208
void libmesh_ignore(const T &)
const std::vector< RealGradient > & get_Sobolev_dweight() const
Definition: fe_base.h:428
const std::vector< std::vector< OutputGradient > > & get_dphi() const
Definition: fe_base.h:216
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
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.
Parameters parameters
Data structure holding arbitrary parameters.
const MeshBase & get_mesh() const
const T_sys & get_system(const std::string &name) const
virtual void add(const numeric_index_type i, const T value)=0
Adds value to each entry of the vector.
This class forms the foundation from which generic finite elements may be derived.
virtual void reinit(const Elem *elem, const std::vector< Point > *const pts=libmesh_nullptr, const std::vector< Real > *const weights=libmesh_nullptr)=0
This is at the core of this class.
const std::vector< OutputGradient > & get_dphase() const
Definition: fe_base.h:404
int main ( int  argc,
char **  argv 
)

Definition at line 83 of file miscellaneous_ex1.C.

References libMesh::EquationSystems::add_system(), assemble_wave(), libMesh::MeshTools::Generation::build_cube(), libMesh::InfElemBuilder::build_inf_elem(), libMesh::LibMeshInit::comm(), libMesh::MeshBase::element_ptr_range(), libMesh::MeshBase::find_neighbors(), libMesh::FIRST, libMesh::EquationSystems::get_system(), libMesh::HEX8, libMesh::TriangleWrapper::init(), libMesh::EquationSystems::init(), mesh, libMesh::out, libMesh::EquationSystems::parameters, libMesh::EquationSystems::print_info(), libMesh::MeshBase::print_info(), libMesh::Real, libMesh::Parameters::set(), libMesh::WRITE, libMesh::ExodusII_IO::write(), and libMesh::EquationSystems::write().

84 {
85  // Initialize libMesh, like in example 2.
86  LibMeshInit init (argc, argv);
87 
88  // This example requires Infinite Elements
89 #ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
90  libmesh_example_requires(false, "--enable-ifem");
91 #else
92 
93  // Skip this 3D example if libMesh was compiled as 1D/2D-only.
94  libmesh_example_requires(3 <= LIBMESH_DIM, "3D support");
95 
96  // Tell the user what we are doing.
97  libMesh::out << "Running ex6 with dim = 3" << std::endl << std::endl;
98 
99  // Create a serialized mesh, distributed across the default MPI
100  // communicator.
101  // InfElemBuilder still requires some updates to be DistributedMesh
102  // compatible
103 
104  ReplicatedMesh mesh(init.comm());
105 
106  // Use the internal mesh generator to create elements
107  // on the square [-1,1]^3, of type Hex8.
109  4, 4, 4,
110  -1., 1.,
111  -1., 1.,
112  -1., 1.,
113  HEX8);
114 
115  // Print information about the mesh to the screen.
116  mesh.print_info();
117 
118  // Write the mesh before the infinite elements are added
119 #ifdef LIBMESH_HAVE_EXODUS_API
120  ExodusII_IO(mesh).write ("orig_mesh.e");
121 #endif
122 
123  // Normally, when a mesh is imported or created in
124  // libMesh, only conventional elements exist. The infinite
125  // elements used here, however, require prescribed
126  // nodal locations (with specified distances from an imaginary
127  // origin) and configurations that a conventional mesh creator
128  // in general does not offer. Therefore, an efficient method
129  // for building infinite elements is offered. It can account
130  // for symmetry planes and creates infinite elements in a fully
131  // automatic way.
132  //
133  // Right now, the simplified interface is used, automatically
134  // determining the origin. Check MeshBase for a generalized
135  // method that can even return the element faces of interior
136  // vibrating surfaces. The bool determines whether to be
137  // verbose.
138  InfElemBuilder builder(mesh);
139  builder.build_inf_elem(true);
140 
141  // Reassign subdomain_id() of all infinite elements.
142  // Otherwise, the exodus-api will fail on the mesh.
143  for (auto & elem : mesh.element_ptr_range())
144  if (elem->infinite())
145  elem->subdomain_id() = 1;
146 
147  // Print information about the mesh to the screen.
148  mesh.print_info();
149 
150  // Write the mesh with the infinite elements added.
151  // Compare this to the original mesh.
152 #ifdef LIBMESH_HAVE_EXODUS_API
153  ExodusII_IO(mesh).write ("ifems_added.e");
154 #endif
155 
156  // After building infinite elements, we have to let
157  // the elements find their neighbors again.
159 
160  // Create an equation systems object, where ThinSystem
161  // offers only the crucial functionality for solving a
162  // system. Use ThinSystem when you want the sleekest
163  // system possible.
164  EquationSystems equation_systems (mesh);
165 
166  // Declare the system and its variables.
167  // Create a system named "Wave". This can
168  // be a simple, steady system
169  equation_systems.add_system<LinearImplicitSystem> ("Wave");
170 
171  // Create an FEType describing the approximation
172  // characteristics of the InfFE object. Note that
173  // the constructor automatically defaults to some
174  // sensible values. But use FIRST order
175  // approximation.
176  FEType fe_type(FIRST);
177 
178  // Add the variable "p" to "Wave". Note that there exist
179  // various approaches in adding variables. In example 3,
180  // add_variable took the order of approximation and used
181  // default values for the FEFamily, while here the FEType
182  // is used.
183  equation_systems.get_system("Wave").add_variable("p", fe_type);
184 
185  // Give the system a pointer to the matrix assembly
186  // function.
187  equation_systems.get_system("Wave").attach_assemble_function (assemble_wave);
188 
189  // Set the speed of sound and fluid density
190  // as EquationSystems parameter,
191  // so that assemble_wave() can access it.
192  equation_systems.parameters.set<Real>("speed") = 1.;
193  equation_systems.parameters.set<Real>("fluid density") = 1.;
194 
195  // Initialize the data structures for the equation system.
196  equation_systems.init();
197 
198  // Prints information about the system to the screen.
199  equation_systems.print_info();
200 
201  // Solve the system "Wave".
202  equation_systems.get_system("Wave").solve();
203 
204  // Write the whole EquationSystems object to file.
205  // For infinite elements, the concept of nodal_soln()
206  // is not applicable. Therefore, writing the mesh in
207  // some format @e always gives all-zero results at
208  // the nodes of the infinite elements. Instead,
209  // use the FEInterface::compute_data() methods to
210  // determine physically correct results within an
211  // infinite element.
212  equation_systems.write ("eqn_sys.dat", WRITE);
213 
214  // All done.
215  return 0;
216 
217 #endif // else part of ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
218 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
This is the EquationSystems class.
The ReplicatedMesh class is derived from the MeshBase class, and is used to store identical copies of...
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.
MeshBase & mesh
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
Definition: libmesh.h:62
virtual void find_neighbors(const bool reset_remote_elements=false, const bool reset_current_list=true)=0
Locate element face (edge in 2D) neighbors.
virtual SimpleRange< element_iterator > element_ptr_range()=0
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
This class is used to build infinite elements on top of an existing mesh.
void assemble_wave(EquationSystems &es, const std::string &system_name)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual void write(const std::string &fname) libmesh_override
This method implements writing a mesh to a specified file.
Definition: exodusII_io.C:789
OStreamProxy out
void print_info(std::ostream &os=libMesh::out) const
Prints relevant information about the mesh.
Definition: mesh_base.C:448
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.