libMesh
Classes | Namespaces | Functions
vector_fe_ex5.C File Reference

Go to the source code of this file.

Classes

class  libMesh::NumericVector< T >
 Provides a uniform interface to vector storage schemes for different linear algebra libraries. More...
 
class  libMesh::SparseMatrix< T >
 Generic sparse matrix. More...
 

Namespaces

 libMesh
 The libMesh namespace provides an interface to certain functionality in the library.
 

Functions

void compute_residual (const NumericVector< Number > &X, NumericVector< Number > &R, NonlinearImplicitSystem &S)
 
void compute_jacobian (const NumericVector< Number > &X, SparseMatrix< Number > &J, NonlinearImplicitSystem &S)
 
int main (int argc, char **argv)
 

Function Documentation

◆ compute_jacobian()

void compute_jacobian ( const NumericVector< Number > &  X,
SparseMatrix< Number > &  J,
NonlinearImplicitSystem S 
)

Definition at line 315 of file assembly.C.

References libMesh::Elem::active(), libMesh::SparseMatrix< T >::add_matrix(), libMesh::FEGenericBase< OutputType >::build(), libMesh::FEType::default_quadrature_order(), dim, libMesh::DofMap::dof_indices(), libMesh::Parameters::get(), libMesh::System::get_dof_map(), libMesh::System::get_equation_systems(), libMesh::System::get_mesh(), libMesh::DofObject::id(), libMesh::FEMap::inverse_map(), libMesh::Elem::level(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::System::name(), libMesh::Elem::neighbor_ptr(), libMesh::EquationSystems::parameters, libMesh::Utility::pow(), libMesh::Real, libMesh::DenseMatrix< T >::resize(), and libMesh::DofMap::variable_type().

Referenced by libMesh::FirstOrderUnsteadySolver::compute_second_order_eqns(), SigmaPhysics::element_time_derivative(), LaplaceSystem::element_time_derivative(), PoissonSystem::element_time_derivative(), HeatSystem::element_time_derivative(), main(), and LaplaceSystem::side_constraint().

318 {
319  // It is a good idea to make sure we are assembling
320  // the proper system.
321  libmesh_assert_equal_to(system.name(), "Poisson");
322 
323  // Get the DG parameters
324  auto & es = system.get_equation_systems();
325  const auto sigma = es.parameters.get<Real>("sigma");
326  const auto epsilon = es.parameters.get<Real>("epsilon");
327 
328  // Get a constant reference to the mesh object.
329  const MeshBase & mesh = system.get_mesh();
330 
331  // The dimension that we are running
332  const auto dim = mesh.mesh_dimension();
333 
334  // A reference to the DofMap object for this system. The DofMap
335  // object handles the index translation from node and element numbers
336  // to degree of freedom numbers. We will talk more about the DofMap
337  // in future examples.
338  const DofMap & dof_map = system.get_dof_map();
339 
340  // Get a constant reference to the Finite Element type
341  // for the first (and only) variable in the system.
342  FEType fe_type = dof_map.variable_type(0);
343 
344  // Build a Finite Element object of the specified type.
345  // Note that FEVectorBase is a typedef for the templated FE
346  // class.
347  std::unique_ptr<FEVectorBase> fe(FEVectorBase::build(dim, fe_type));
348 
349  // An automatically determined Gauss quadrature rule for numerical integration.
350  QGauss qrule(dim, fe_type.default_quadrature_order());
351 
352  // Tell the finite element object to use our quadrature rule.
353  fe->attach_quadrature_rule(&qrule);
354 
355  // Declare a special finite element object for face integration.
356  std::unique_ptr<FEVectorBase> fe_face(FEVectorBase::build(dim, fe_type));
357  // And for neighbor integration
358  std::unique_ptr<FEVectorBase> fe_neighbor_face(FEVectorBase::build(dim, fe_type));
359 
360  // Boundary integration requires one quadrature rule,
361  // with dimensionality one less than the dimensionality
362  // of the element.
363  QGauss qface(dim - 1, fe_type.default_quadrature_order());
364 
365  // Tell the face finite element objects to use our
366  // quadrature rule.
367  fe_face->attach_quadrature_rule(&qface);
368  fe_neighbor_face->attach_quadrature_rule(&qface);
369 
370  // Here we define some references to cell-specific data that
371  // will be used to assemble the linear system.
372  //
373  // The element Jacobian * quadrature weight at each integration point.
374  const auto & JxW = fe->get_JxW();
375 
376  // The element shape function gradients evaluated at the quadrature points.
377  const auto & dphi = fe->get_dphi();
378 
379  // face integration points
380  const auto & face_xyz = fe_face->get_xyz();
381 
382  // Face shape function values
383  const auto & phi_face = fe_face->get_phi();
384 
385  // Face shape function gradients
386  const auto & dphi_face = fe_face->get_dphi();
387 
388  // Neighbor shape function values
389  const auto & phi_neighbor = fe_neighbor_face->get_phi();
390 
391  // Neighbor face shape function gradients
392  const auto & dphi_neighbor = fe_neighbor_face->get_dphi();
393 
394  // face normals
395  const auto & normals = fe_face->get_normals();
396 
397  // face JxW
398  const auto & JxW_face = fe_face->get_JxW();
399 
404 
405  // This vector will hold the degree of freedom indices for
406  // the element. These define where in the global system
407  // the element degrees of freedom get mapped.
408  std::vector<dof_id_type> dof_indices;
409  std::vector<dof_id_type> dof_indices_neighbor;
410 
411  // To avoid extraneous allocation when building element sides (to compute their volumes)
412  ElemSideBuilder side_builder;
413 
414  // Now we will loop over all the elements in the mesh.
415  // We will compute the element matrix and right-hand-side
416  // contribution.
417  //
418  // Element iterators are a nice way to iterate through all the
419  // elements, or all the elements that have some property. The
420  // iterator el will iterate from the first to the last element on
421  // the local processor. The iterator end_el tells us when to stop.
422  // It is smart to make this one const so that we don't accidentally
423  // mess it up! In case users later modify this program to include
424  // refinement, we will be safe and will only consider the active
425  // elements; hence we use a variant of the active_elem_iterator.
426  for (const auto & elem : mesh.active_local_element_ptr_range())
427  {
428  // Get the degree of freedom indices for the
429  // current element. These define where in the global
430  // matrix and right-hand-side this element will
431  // contribute to.
432  dof_map.dof_indices(elem, dof_indices);
433 
434  // Compute the element-specific data for the current
435  // element. This involves computing the location of the
436  // quadrature points (q_point) and the shape functions
437  // (phi, dphi) for the current element.
438  fe->reinit(elem);
439 
440  libmesh_assert_msg(dphi.size() == dof_indices.size(),
441  "dphi size doesn't match dof_indices size");
442 
443  // DenseMatrix::resize() member will automatically zero out the matrix.
444  Kee.resize(dof_indices.size(), dof_indices.size());
445 
446  // diffusion elemental jacobian
447  for (unsigned int qp = 0; qp < qrule.n_points(); qp++)
448  for (std::size_t i = 0; i < dof_indices.size(); i++)
449  for (std::size_t j = 0; j < dof_indices.size(); j++)
450  Kee(i, j) += JxW[qp] * dphi[i][qp].contract(dphi[j][qp]);
451 
452  // Now we consider jacobian contributions from the sides
453  for (auto side : elem->side_index_range())
454  {
455  // We need to compute h for penalty terms
456  const auto side_volume = side_builder(*elem, side).volume();
457  fe_face->reinit(elem, side);
458  const auto elem_b_order = static_cast<unsigned int>(fe_face->get_order());
459  const auto h_elem = elem->volume() / side_volume / std::pow(elem_b_order, 2);
460 
461  // No neighbor means we must be on a boundary
462  if (!elem->neighbor_ptr(side))
463  {
464  for (unsigned int qp = 0; qp < qface.n_points(); qp++)
465  for (std::size_t i = 0; i < dof_indices.size(); i++)
466  for (std::size_t j = 0; j < dof_indices.size(); j++)
467 
468  {
469  Kee(i, j) -= dphi_face[j][qp] * normals[qp] * phi_face[i][qp] * JxW_face[qp];
470  Kee(i, j) +=
471  epsilon * phi_face[j][qp] * dphi_face[i][qp] * normals[qp] * JxW_face[qp];
472  Kee(i, j) += sigma / h_elem * phi_face[j][qp] * phi_face[i][qp] * JxW_face[qp];
473  }
474  }
475  else // We must be on an interior side
476  {
477  const Elem * neighbor = elem->neighbor_ptr(side);
478 
479  const auto elem_id = elem->id();
480  const auto neighbor_id = neighbor->id();
481 
482  // We don't want to erroneously add multiple contributions from the same interior face
483  if ((neighbor->active() && (neighbor->level() == elem->level()) &&
484  (elem_id < neighbor_id)) ||
485  (neighbor->level() < elem->level()))
486  {
487  dof_map.dof_indices(neighbor, dof_indices_neighbor);
488 
489  // Make sure we have the matching quadrature points on face and neighbor
490  std::vector<Point> neighbor_xyz;
491  FEMap::inverse_map(elem->dim(), neighbor, face_xyz,
492  neighbor_xyz);
493 
494  fe_neighbor_face->reinit(neighbor, &neighbor_xyz);
495 
496  libmesh_assert_msg(dphi_neighbor.size() == dof_indices_neighbor.size(),
497  "dphi_neighbor size doesn't match dof_indices_neighbor size");
498 
499  Ken.resize(dof_indices.size(), dof_indices_neighbor.size());
500  Kne.resize(dof_indices_neighbor.size(), dof_indices.size());
501  Knn.resize(dof_indices_neighbor.size(), dof_indices_neighbor.size());
502 
503  // Now add the DG contribution to the local jacobian
504  for (unsigned int qp = 0; qp < qface.n_points(); qp++)
505  {
506  // element-element contribution
507  for (std::size_t i = 0; i < dof_indices.size(); i++)
508  for (std::size_t j = 0; j < dof_indices.size(); j++)
509  {
510  Kee(i, j) -= 0.5 * dphi_face[j][qp] * normals[qp] * phi_face[i][qp] * JxW_face[qp];
511  Kee(i, j) +=
512  epsilon * 0.5 * phi_face[j][qp] * dphi_face[i][qp] * normals[qp] * JxW_face[qp];
513  Kee(i, j) += sigma / h_elem * phi_face[j][qp] * phi_face[i][qp] * JxW_face[qp];
514  }
515  // element-neighbor contribution
516  for (std::size_t i = 0; i < dof_indices.size(); i++)
517  for (std::size_t j = 0; j < dof_indices_neighbor.size(); j++)
518  {
519  Ken(i, j) -=
520  0.5 * dphi_neighbor[j][qp] * normals[qp] * phi_face[i][qp] * JxW_face[qp];
521  Ken(i, j) += epsilon * 0.5 * -phi_neighbor[j][qp] * dphi_face[i][qp] * normals[qp] *
522  JxW_face[qp];
523  Ken(i, j) += sigma / h_elem * -phi_neighbor[j][qp] * phi_face[i][qp] * JxW_face[qp];
524  }
525  // Neighbor-element contribution
526  for (std::size_t i = 0; i < dof_indices_neighbor.size(); i++)
527  for (std::size_t j = 0; j < dof_indices_neighbor.size(); j++)
528  {
529  Kne(i, j) +=
530  0.5 * dphi_face[j][qp] * normals[qp] * phi_neighbor[i][qp] * JxW_face[qp];
531  Kne(i, j) -= epsilon * 0.5 * phi_face[j][qp] * dphi_neighbor[i][qp] * normals[qp] *
532  JxW_face[qp];
533  Kne(i, j) -= sigma / h_elem * phi_face[j][qp] * phi_neighbor[i][qp] * JxW_face[qp];
534  }
535  // Neighbor-neighbor contribution
536  for (std::size_t i = 0; i < dof_indices_neighbor.size(); i++)
537  for (std::size_t j = 0; j < dof_indices_neighbor.size(); j++)
538  {
539  Knn(i, j) +=
540  0.5 * dphi_neighbor[j][qp] * normals[qp] * phi_neighbor[i][qp] * JxW_face[qp];
541  Knn(i, j) -= epsilon * 0.5 * -phi_neighbor[j][qp] * dphi_neighbor[i][qp] *
542  normals[qp] * JxW_face[qp];
543  Knn(i, j) -=
544  sigma / h_elem * -phi_neighbor[j][qp] * phi_neighbor[i][qp] * JxW_face[qp];
545  }
546  }
547 
548  J.add_matrix(Ken, dof_indices, dof_indices_neighbor);
549  J.add_matrix(Kne, dof_indices_neighbor, dof_indices);
550  J.add_matrix(Knn, dof_indices_neighbor, dof_indices_neighbor);
551  } // whether we've done this internal side integral before
552  } // whether we're on an internal side
553  } // loop over sides
554 
555  J.add_matrix(Kee, dof_indices, dof_indices);
556  } // element loop
557 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:196
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Definition: dof_map.C:2164
unsigned int dim
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:2220
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
MeshBase & mesh
Order default_quadrature_order() const
Definition: fe_type.h:371
This is the MeshBase class.
Definition: mesh_base.h:75
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:179
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.
T pow(const T &x)
Definition: utility.h:328
dof_id_type id() const
Definition: dof_object.h:828
Helper for building element sides that minimizes the construction of new elements.
const Elem * neighbor_ptr(unsigned int i) const
Definition: elem.h:2598
unsigned int level() const
Definition: elem.h:3074
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void resize(const unsigned int new_m, const unsigned int new_n)
Resizes the matrix to the specified size and calls zero().
Definition: dense_matrix.h:895
This class implements specific orders of Gauss quadrature.
unsigned int mesh_dimension() const
Definition: mesh_base.C:372
bool active() const
Definition: elem.h:2941

◆ compute_residual()

void compute_residual ( const NumericVector< Number > &  X,
NumericVector< Number > &  R,
NonlinearImplicitSystem S 
)

Vectors to hold the local solution degree of freedom values

Vector to hold the local solution

Definition at line 27 of file assembly.C.

References libMesh::Elem::active(), libMesh::NumericVector< T >::add_vector(), libMesh::FEGenericBase< OutputType >::build(), libMesh::FEType::default_quadrature_order(), dim, libMesh::DofMap::dof_indices(), exact_solution(), forcing_function(), libMesh::Parameters::get(), libMesh::NumericVector< T >::get(), libMesh::System::get_dof_map(), libMesh::System::get_equation_systems(), libMesh::System::get_mesh(), libMesh::DofObject::id(), libMesh::FEMap::inverse_map(), libMesh::Elem::level(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::System::name(), libMesh::Elem::neighbor_ptr(), libMesh::EquationSystems::parameters, libMesh::Utility::pow(), libMesh::Real, libMesh::DenseVector< T >::resize(), and libMesh::DofMap::variable_type().

Referenced by main().

30 {
31  // It is a good idea to make sure we are assembling
32  // the proper system.
33  libmesh_assert_equal_to(system.name(), "Poisson");
34 
35  // Get the DG parameters
36  auto & es = system.get_equation_systems();
37  const auto sigma = es.parameters.get<Real>("sigma");
38  const auto epsilon = es.parameters.get<Real>("epsilon");
39 
40  // Get a constant reference to the mesh object.
41  const MeshBase & mesh = system.get_mesh();
42 
43  // The dimension that we are running
44  const auto dim = mesh.mesh_dimension();
45 
46  // A reference to the DofMap object for this system. The DofMap
47  // object handles the index translation from node and element numbers
48  // to degree of freedom numbers. We will talk more about the DofMap
49  // in future examples.
50  const DofMap & dof_map = system.get_dof_map();
51 
52  // Get a constant reference to the Finite Element type
53  // for the first (and only) variable in the system.
54  FEType fe_type = dof_map.variable_type(0);
55 
56  // Build a Finite Element object of the specified type.
57  // Note that FEVectorBase is a typedef for the templated FE
58  // class.
59  std::unique_ptr<FEVectorBase> fe(FEVectorBase::build(dim, fe_type));
60 
61  // An automatically determined Gauss quadrature rule for numerical integration.
62  QGauss qrule(dim, fe_type.default_quadrature_order());
63 
64  // Tell the finite element object to use our quadrature rule.
65  fe->attach_quadrature_rule(&qrule);
66 
67  // Declare a special finite element object for face integration.
68  std::unique_ptr<FEVectorBase> fe_face(FEVectorBase::build(dim, fe_type));
69  // And for neighbor integration
70  std::unique_ptr<FEVectorBase> fe_neighbor_face(FEVectorBase::build(dim, fe_type));
71 
72  // Boundary integration requires one quadrature rule,
73  // with dimensionality one less than the dimensionality
74  // of the element.
75  QGauss qface(dim - 1, fe_type.default_quadrature_order());
76 
77  // Tell the face finite element objects to use our
78  // quadrature rule.
79  fe_face->attach_quadrature_rule(&qface);
80  fe_neighbor_face->attach_quadrature_rule(&qface);
81 
82  // Here we define some references to cell-specific data that
83  // will be used to assemble the linear system.
84  //
85  // The element Jacobian * quadrature weight at each integration point.
86  const auto & JxW = fe->get_JxW();
87 
88  // element integration points
89  const auto & xyz = fe->get_xyz();
90 
91  // The element shape function values evaluated at the quadrature points
92  const auto & phi = fe->get_phi();
93 
94  // The element shape function gradients evaluated at the quadrature points.
95  const auto & dphi = fe->get_dphi();
96 
97  // face integration points
98  const auto & face_xyz = fe_face->get_xyz();
99 
100  // Face shape function values
101  const auto & phi_face = fe_face->get_phi();
102 
103  // Face shape function gradients
104  const auto & dphi_face = fe_face->get_dphi();
105 
106  // Neighbor shape function values
107  const auto & phi_neighbor = fe_neighbor_face->get_phi();
108 
109  // Neighbor face shape function gradients
110  const auto & dphi_neighbor = fe_neighbor_face->get_dphi();
111 
112  // face normals
113  const auto & normals = fe_face->get_normals();
114 
115  // face JxW
116  const auto & JxW_face = fe_face->get_JxW();
117 
120 
121  // This vector will hold the degree of freedom indices for
122  // the element. These define where in the global system
123  // the element degrees of freedom get mapped.
124  std::vector<dof_id_type> dof_indices;
125  std::vector<dof_id_type> dof_indices_neighbor;
126 
128  std::vector<Number> dof_u;
129  std::vector<Number> dof_u_neighbor;
130 
132  std::vector<VectorValue<Number>> u;
133  std::vector<TensorValue<Number>> grad_u;
134  std::vector<VectorValue<Number>> u_neighbor;
135  std::vector<TensorValue<Number>> grad_u_neighbor;
136 
137  // Now we will loop over all the elements in the mesh.
138  // We will compute the element matrix and right-hand-side
139  // contribution.
140  //
141  // Element iterators are a nice way to iterate through all the
142  // elements, or all the elements that have some property. The
143  // iterator el will iterate from the first to the last element on
144  // the local processor. The iterator end_el tells us when to stop.
145  // It is smart to make this one const so that we don't accidentally
146  // mess it up! In case users later modify this program to include
147  // refinement, we will be safe and will only consider the active
148  // elements; hence we use a variant of the active_elem_iterator.
149  for (const auto & elem : mesh.active_local_element_ptr_range())
150  {
151  // Get the degree of freedom indices for the
152  // current element. These define where in the global
153  // matrix and right-hand-side this element will
154  // contribute to.
155  dof_map.dof_indices(elem, dof_indices);
156 
157  // Compute the element-specific data for the current
158  // element. This involves computing the location of the
159  // quadrature points (q_point) and the shape functions
160  // (phi, dphi) for the current element.
161  fe->reinit(elem);
162 
163  libmesh_assert_msg(dphi.size() == dof_indices.size(),
164  "dphi size doesn't match dof_indices size");
165 
166  // DenseVector::resize() member will automatically zero out the vector.
167  Fe.resize(dof_indices.size());
168 
169  // Get the local solution vector
170  X.get(dof_indices, dof_u);
171 
172  // build the element solution gradient
173  grad_u.resize(qrule.n_points());
174  for (unsigned int qp = 0; qp < qrule.n_points(); qp++)
175  {
176  grad_u[qp] = 0;
177  for (std::size_t i = 0; i < dof_indices.size(); i++)
178  grad_u[qp] += dof_u[i] * dphi[i][qp];
179  }
180 
181  for (unsigned int qp = 0; qp < qrule.n_points(); qp++)
182  {
183  auto forcing_func = VectorValue<Number>(forcing_function(0, xyz[qp](0), xyz[qp](1)),
184  forcing_function(1, xyz[qp](0), xyz[qp](1)),
185  0);
186  for (std::size_t i = 0; i < dof_indices.size(); i++)
187  {
188  // diffusion elemental residual
189  Fe(i) += JxW[qp] * dphi[i][qp].contract(grad_u[qp]);
190 
191  // forcing function
192  Fe(i) += JxW[qp] * phi[i][qp] * forcing_func;
193  }
194  }
195 
196  // To avoid extraneous allocation when building element sides (to compute their volumes)
197  ElemSideBuilder side_builder;
198 
199  // Now we consider residual contributions from the sides
200  for (auto side : elem->side_index_range())
201  {
202  // We need to compute h for penalty terms
203  const auto side_volume = side_builder(*elem, side).volume();
204  fe_face->reinit(elem, side);
205  const auto elem_b_order = static_cast<unsigned int>(fe_face->get_order());
206  const auto h_elem = elem->volume() / side_volume / std::pow(elem_b_order, 2);
207 
208  // build the face solution value and gradient
209  u.resize(qface.n_points());
210  grad_u.resize(qface.n_points());
211  for (unsigned int qp = 0; qp < qface.n_points(); qp++)
212  {
213  u[qp] = 0;
214  grad_u[qp] = 0;
215  for (std::size_t i = 0; i < dof_indices.size(); i++)
216  {
217  u[qp] += dof_u[i] * phi_face[i][qp];
218  grad_u[qp] += dof_u[i] * dphi_face[i][qp];
219  }
220  }
221 
222  // No neighbor means we must be on a boundary
223  if (!elem->neighbor_ptr(side))
224  {
225  for (unsigned int qp = 0; qp < qface.n_points(); qp++)
226  {
227  auto fn = VectorValue<Number>(exact_solution(0, face_xyz[qp](0), face_xyz[qp](1)),
228  exact_solution(1, face_xyz[qp](0), face_xyz[qp](1)),
229  0);
230  for (std::size_t i = 0; i < dof_indices.size(); i++)
231  {
232  Fe(i) -= grad_u[qp] * normals[qp] * phi_face[i][qp] * JxW_face[qp];
233  Fe(i) += epsilon * (u[qp] - fn) * dphi_face[i][qp] * normals[qp] * JxW_face[qp];
234  Fe(i) += sigma / h_elem * (u[qp] - fn) * phi_face[i][qp] * JxW_face[qp];
235  }
236  }
237  }
238  else // We must be on an interior side
239  {
240  const Elem * neighbor = elem->neighbor_ptr(side);
241 
242  const auto elem_id = elem->id();
243  const auto neighbor_id = neighbor->id();
244 
245  // We don't want to erroneously add multiple contributions from the same interior face
246  if ((neighbor->active() && (neighbor->level() == elem->level()) &&
247  (elem_id < neighbor_id)) ||
248  (neighbor->level() < elem->level()))
249  {
250  dof_map.dof_indices(neighbor, dof_indices_neighbor);
251 
252  // Make sure we have the matching quadrature points on face and neighbor
253  std::vector<Point> neighbor_xyz;
254  FEMap::inverse_map(elem->dim(), neighbor, face_xyz,
255  neighbor_xyz);
256 
257  fe_neighbor_face->reinit(neighbor, &neighbor_xyz);
258 
259  libmesh_assert_msg(dphi_neighbor.size() == dof_indices_neighbor.size(),
260  "dphi_neighbor size doesn't match dof_indices_neighbor size");
261 
262  Fn.resize(dof_indices_neighbor.size());
263  X.get(dof_indices_neighbor, dof_u_neighbor);
264 
265  // build the neighbor solution value and gradient
266  u_neighbor.resize(qface.n_points());
267  grad_u_neighbor.resize(qface.n_points());
268  for (unsigned int qp = 0; qp < qface.n_points(); qp++)
269  {
270  u_neighbor[qp] = 0;
271  grad_u_neighbor[qp] = 0;
272  for (std::size_t i = 0; i < dof_indices_neighbor.size(); i++)
273  {
274  u_neighbor[qp] += dof_u_neighbor[i] * phi_neighbor[i][qp];
275  grad_u_neighbor[qp] += dof_u_neighbor[i] * dphi_neighbor[i][qp];
276  }
277  }
278 
279  // Now add the DG contribution to the local residual
280  for (unsigned int qp = 0; qp < qface.n_points(); qp++)
281  {
282  // element contribution
283  for (std::size_t i = 0; i < dof_indices.size(); i++)
284  {
285  Fe(i) -= 0.5 * (grad_u[qp] * normals[qp] + grad_u_neighbor[qp] * normals[qp]) *
286  phi_face[i][qp] * JxW_face[qp];
287  Fe(i) += epsilon * 0.5 * (u[qp] - u_neighbor[qp]) * dphi_face[i][qp] * normals[qp] *
288  JxW_face[qp];
289  Fe(i) += sigma / h_elem * (u[qp] - u_neighbor[qp]) * phi_face[i][qp] * JxW_face[qp];
290  }
291  // Neighbor contribution
292  for (std::size_t i = 0; i < dof_indices_neighbor.size(); i++)
293  {
294  Fn(i) += 0.5 * (grad_u[qp] * normals[qp] + grad_u_neighbor[qp] * normals[qp]) *
295  phi_neighbor[i][qp] * JxW_face[qp];
296  Fn(i) -= epsilon * 0.5 * (u[qp] - u_neighbor[qp]) * dphi_neighbor[i][qp] *
297  normals[qp] * JxW_face[qp];
298  Fn(i) -=
299  sigma / h_elem * (u[qp] - u_neighbor[qp]) * phi_neighbor[i][qp] * JxW_face[qp];
300  }
301  }
302 
303  R.add_vector(Fn, dof_indices_neighbor);
304  } // whether we've done this internal side integral before
305  } // whether we're on an internal side
306  } // loop over sides
307 
308  R.add_vector(Fe, dof_indices);
309  } // element loop
310 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:196
Real exact_solution(const int component, const Real x, const Real y, const Real z=0.)
This is the exact solution that we are trying to obtain.
virtual void get(const std::vector< numeric_index_type > &index, T *values) const
Access multiple components at once.
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Definition: dof_map.C:2164
unsigned int dim
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:396
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]...
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:2220
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
MeshBase & mesh
Order default_quadrature_order() const
Definition: fe_type.h:371
Real forcing_function(const int component, const Real x, const Real y, const Real z=0.)
This class defines a vector in LIBMESH_DIM dimensional Real or Complex space.
This is the MeshBase class.
Definition: mesh_base.h:75
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:179
T pow(const T &x)
Definition: utility.h:328
dof_id_type id() const
Definition: dof_object.h:828
Helper for building element sides that minimizes the construction of new elements.
const Elem * neighbor_ptr(unsigned int i) const
Definition: elem.h:2598
unsigned int level() const
Definition: elem.h:3074
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
This class implements specific orders of Gauss quadrature.
unsigned int mesh_dimension() const
Definition: mesh_base.C:372
bool active() const
Definition: elem.h:2941

◆ main()

int main ( int  argc,
char **  argv 
)

Definition at line 63 of file vector_fe_ex5.C.

References libMesh::EquationSystems::add_system(), libMesh::System::add_variable(), libMesh::MeshTools::Generation::build_square(), libMesh::EquationSystems::build_variable_names(), compute_jacobian(), compute_residual(), libMesh::default_solver_package(), libMesh::TriangleWrapper::init(), libMesh::EquationSystems::init(), mesh, libMesh::MONOMIAL_VEC, libMesh::out, libMesh::EquationSystems::parameters, libMesh::PETSC_SOLVERS, libMesh::EquationSystems::print_info(), libMesh::MeshBase::print_info(), libMesh::QUAD4, libMesh::Real, libMesh::Parameters::set(), libMesh::ExodusII_IO::set_output_variables(), libMesh::TRILINOS_SOLVERS, libMesh::ExodusII_IO::write(), libMesh::ExodusII_IO::write_element_data(), and libMesh::ExodusII_IO::write_equation_systems().

64 {
65  // Initialize libraries.
66  LibMeshInit init(argc, argv);
67 
68  // This example requires a non-linear solver package.
69  libmesh_example_requires(libMesh::default_solver_package() == PETSC_SOLVERS ||
71  "--enable-petsc or --enable-trilinos");
72 
73  // Parse input file
74  GetPot input_file("vector_fe_ex5.in");
75 
76  // But allow the command line to override it.
77  input_file.parse_command_line(argc, argv);
78 
79  // Read DG parameters
80  const Real epsilon = input_file("epsilon", -1);
81  const Real sigma = input_file("sigma", 6);
82 
83  // Read mesh size
84  const std::size_t nx = input_file("nx", 15);
85  const std::size_t ny = input_file("ny", 15);
86 
87  // Read approximation order - default to FIRST
88  const unsigned int approx_order = input_file("approx_order", 1);
89 
90  // Brief message to the user regarding the program name
91  // and command line arguments.
92  libMesh::out << "Running " << argv[0];
93 
94  for (int i = 1; i < argc; i++)
95  libMesh::out << " " << argv[i];
96 
97  libMesh::out << std::endl << std::endl;
98 
99  // Skip this 2D example if libMesh was compiled as 1D-only.
100  libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");
101 
102  // Create a mesh, with dimension to be overridden later, on the
103  // default MPI communicator.
104  Mesh mesh(init.comm());
105 
106  // Use the MeshTools::Generation mesh generator to create a uniform
107  // 2D grid on the square [-1,1]^2. We instruct the mesh generator
108  // to build a mesh of 15x15 QUAD4 elements. Note that at the end of this
109  // function call, the mesh will be prepared and remote elements will be deleted
110  MeshTools::Generation::build_square(mesh, nx, ny, -1., 1., -1., 1., QUAD4);
111 
112  // Print information about the mesh to the screen.
113  mesh.print_info();
114 
115  // Create an equation systems object.
116  EquationSystems equation_systems(mesh);
117  equation_systems.parameters.set<Real>("epsilon") = epsilon;
118  equation_systems.parameters.set<Real>("sigma") = sigma;
119 
120  // Declare the Poisson system and its variables.
121  // The Poisson system is another example of a steady system.
122  auto & nl_system = equation_systems.add_system<NonlinearImplicitSystem>("Poisson");
123 
124  // Adds the variable "u" to "Poisson". Variable "u" will be approximated using
125  // first-order (can be overridden with command line) vector monomial elements.
126  // Since the mesh is 2-D, "u" will have two components.
127  nl_system.add_variable("u", static_cast<Order>(approx_order), MONOMIAL_VEC);
128 
129  // Set the residual and Jacobian evaluation functions
130  auto & nl_solver = *nl_system.nonlinear_solver;
131  nl_solver.residual = compute_residual;
132  nl_solver.jacobian = compute_jacobian;
133 
134  // Initialize the data structures for the equation system.
135  equation_systems.init();
136 
137  // Prints information about the system to the screen.
138  equation_systems.print_info();
139 
140  // A p=0 "solve" here is a hack to test ExodusII output; the
141  // solve() is stymied by division by zero.
142  std::unique_ptr<FPEDisabler> maybe_disable_fpe;
143  if (! approx_order) /* CONSTANT, MONOMIAL_VEC */
144  maybe_disable_fpe = std::make_unique<FPEDisabler>();
145 
146  nl_system.solve();
147 
148 #ifdef LIBMESH_HAVE_EXODUS_API
149  if (! approx_order) /* CONSTANT, MONOMIAL_VEC */
150  {
151  // Warn about failed solve for CONSTANT approximations (Poisson needs at least p=1)
152  libMesh::out << "\nWARNING: The elemental vector order is CONSTANT. The solution" << std::endl
153  << "written out to 'out_constant.e' will be trivial." << std::endl;
154 
155  // Need to get string of variable names (vector components) for 'set_output_variables()'
156  std::vector<std::string> varnames;
157  equation_systems.build_variable_names(varnames);
158 
159  // Create ExodusII object with an unambiguous filename (indicating this solution is special)
160  ExodusII_IO exo_ptr(mesh);
161  exo_ptr.write("out_constant.e");
162  exo_ptr.set_output_variables(varnames);
163  exo_ptr.write_element_data(equation_systems);
164  }
165  else
166  ExodusII_IO(mesh).write_equation_systems("out.e", equation_systems);
167 #endif
168 
169  // All done.
170  return 0;
171 }
This is the EquationSystems class.
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
Definition: exodusII_io.h:52
MeshBase & mesh
void build_square(UnstructuredMesh &mesh, const unsigned int nx, const unsigned int ny, const Real xmin=0., const Real xmax=1., const Real ymin=0., const Real ymax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
A specialized build_cube() for 2D meshes.
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
Definition: libmesh.h:90
void compute_jacobian(const NumericVector< Number > &X, SparseMatrix< Number > &J, NonlinearImplicitSystem &S)
Definition: assembly.C:315
void compute_residual(const NumericVector< Number > &X, NumericVector< Number > &R, NonlinearImplicitSystem &S)
Definition: assembly.C:27
SolverPackage default_solver_package()
Definition: libmesh.C:1117
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.
Manages consistently variables, degrees of freedom, coefficient vectors, matrices and non-linear solv...
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.
Definition: system.C:1357
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
OStreamProxy out
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50