libMesh
Classes | Functions
systems_test.C File Reference

Go to the source code of this file.

Classes

class  AugmentSparsityOnNodes
 
struct  TripleFunction
 
class  SystemsTest
 

Functions

void assemble_matrix_and_rhs (EquationSystems &es, const std::string &)
 
void assembly_with_dg_fem_context (EquationSystems &es, const std::string &)
 
Number cubic_test (const Point &p, const Parameters &, const std::string &, const std::string &)
 
Number new_linear_test (const Point &p, const Parameters &, const std::string &, const std::string &)
 
Number disc_thirds_test (const Point &p, const Parameters &, const std::string &, const std::string &)
 
 CPPUNIT_TEST_SUITE_REGISTRATION (SystemsTest)
 

Function Documentation

◆ assemble_matrix_and_rhs()

void assemble_matrix_and_rhs ( EquationSystems es,
const std::string &   
)

Definition at line 123 of file systems_test.C.

References libMesh::SparseMatrix< T >::add_matrix(), libMesh::NumericVector< T >::add_vector(), libMesh::NumericVector< T >::close(), libMesh::SparseMatrix< T >::close(), libMesh::DofObject::dof_number(), libMesh::System::get_dof_map(), libMesh::EquationSystems::get_mesh(), libMesh::EquationSystems::get_system(), libMesh::ImplicitSystem::get_system_matrix(), mesh, libMesh::MeshBase::node_ref(), libMesh::NODEELEM, libMesh::System::number(), libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), libMesh::ExplicitSystem::rhs, libMesh::Elem::type(), and libMesh::System::variable_number().

Referenced by ConstraintOperatorTest::test1DCoarseningNewNodes(), ConstraintOperatorTest::test1DCoarseningOperator(), ConstraintOperatorTest::testCoreform(), and SystemsTest::testDofCouplingWithVarGroups().

125 {
126  const MeshBase& mesh = es.get_mesh();
128  const DofMap& dof_map = system.get_dof_map();
129 
132 
133  std::vector<dof_id_type> dof_indices;
134 
135  SparseMatrix<Number> & matrix = system.get_system_matrix();
136 
137  MeshBase::const_element_iterator el = mesh.active_local_elements_begin();
138  const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end();
139 
140  for ( ; el != end_el; ++el)
141  {
142  const Elem* elem = *el;
143 
144  if(elem->type() == NODEELEM)
145  {
146  continue;
147  }
148 
149  dof_map.dof_indices (elem, dof_indices);
150  const unsigned int n_dofs = dof_indices.size();
151 
152  Ke.resize (n_dofs, n_dofs);
153  Fe.resize (n_dofs);
154 
155  for(unsigned int i=0; i<n_dofs; i++)
156  {
157  Ke(i,i) = 1.;
158  Fe(i) = 1.;
159  }
160 
161  matrix.add_matrix (Ke, dof_indices);
162  system.rhs->add_vector (Fe, dof_indices);
163  }
164 
165  // Add matrix for extra coupled dofs
166  {
167  const Node & node_1 = mesh.node_ref(1);
168  const Node & node_2 = mesh.node_ref(2);
169  dof_indices.resize(6);
170  dof_indices[0] =
171  node_1.dof_number(system.number(), system.variable_number("u"), 0);
172  dof_indices[1] =
173  node_1.dof_number(system.number(), system.variable_number("v"), 0);
174  dof_indices[2] =
175  node_1.dof_number(system.number(), system.variable_number("w"), 0);
176 
177  dof_indices[3] =
178  node_2.dof_number(system.number(), system.variable_number("u"), 0);
179  dof_indices[4] =
180  node_2.dof_number(system.number(), system.variable_number("v"), 0);
181  dof_indices[5] =
182  node_2.dof_number(system.number(), system.variable_number("w"), 0);
183 
184  const unsigned int n_dofs = dof_indices.size();
185  Ke.resize (n_dofs, n_dofs);
186  Fe.resize (n_dofs);
187 
188  for(unsigned int i=0; i<n_dofs; i++)
189  {
190  Ke(i,i) = 1.;
191  Fe(i) = 1.;
192  }
193 
194  matrix.add_matrix (Ke, dof_indices);
195  system.rhs->add_vector (Fe, dof_indices);
196  }
197 
198  system.rhs->close();
199  matrix.close();
200 }
dof_id_type dof_number(const unsigned int s, const unsigned int var, const unsigned int comp) const
Definition: dof_object.h:1032
A Node is like a Point, but with more information.
Definition: node.h:52
The definition of the const_element_iterator struct.
Definition: mesh_base.h:2216
Manages consistently variables, degrees of freedom, coefficient vectors, matrices and linear solvers ...
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]...
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
MeshBase & mesh
NumericVector< Number > * rhs
The system matrix.
const T_sys & get_system(std::string_view name) const
This is the MeshBase class.
Definition: mesh_base.h:75
unsigned int variable_number(std::string_view var) const
Definition: system.C:1609
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:179
unsigned int number() const
Definition: system.h:2350
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.
virtual void close()=0
Calls the NumericVector&#39;s internal assembly routines, ensuring that the values are consistent across ...
virtual void close()=0
Calls the SparseMatrix&#39;s internal assembly routines, ensuring that the values are consistent across p...
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().
Definition: dense_matrix.h:895
virtual const Node & node_ref(const dof_id_type i) const
Definition: mesh_base.h:596
const DofMap & get_dof_map() const
Definition: system.h:2374
virtual ElemType type() const =0
const SparseMatrix< Number > & get_system_matrix() const

◆ assembly_with_dg_fem_context()

void assembly_with_dg_fem_context ( EquationSystems es,
const std::string &   
)

Definition at line 203 of file systems_test.C.

References libMesh::SparseMatrix< T >::add_matrix(), libMesh::NumericVector< T >::add_vector(), libMesh::FEMContext::elem_fe_reinit(), libMesh::DiffContext::get_dof_indices(), libMesh::FEGenericBase< OutputType >::get_dphi(), libMesh::DGFEMContext::get_elem_elem_jacobian(), libMesh::DiffContext::get_elem_jacobian(), libMesh::DGFEMContext::get_elem_neighbor_jacobian(), libMesh::DiffContext::get_elem_residual(), libMesh::FEMContext::get_element_fe(), libMesh::FEMContext::get_element_qrule(), libMesh::FEAbstract::get_JxW(), libMesh::EquationSystems::get_mesh(), libMesh::DGFEMContext::get_neighbor_dof_indices(), libMesh::DGFEMContext::get_neighbor_elem_jacobian(), libMesh::DGFEMContext::get_neighbor_neighbor_jacobian(), libMesh::DGFEMContext::get_neighbor_side_fe(), libMesh::FEGenericBase< OutputType >::get_phi(), libMesh::FEMContext::get_side(), libMesh::FEMContext::get_side_fe(), libMesh::FEMContext::get_side_qrule(), libMesh::EquationSystems::get_system(), libMesh::ImplicitSystem::get_system_matrix(), mesh, libMesh::QBase::n_points(), libMesh::Elem::neighbor_ptr(), libMesh::DGFEMContext::neighbor_side_fe_reinit(), libMesh::FEMContext::pre_fe_reinit(), libMesh::ExplicitSystem::rhs, libMesh::DGFEMContext::set_neighbor(), libMesh::FEMContext::side, and libMesh::DGFEMContext::side_fe_reinit().

Referenced by SystemsTest::testAssemblyWithDgFemContext(), and SystemsTest::testSetSystemParameterOverEquationSystem().

205 {
206  const MeshBase& mesh = es.get_mesh();
208 
211 
212  std::vector<dof_id_type> dof_indices;
213  SparseMatrix<Number> & matrix = system.get_system_matrix();
214 
215  DGFEMContext context(system);
216  {
217  // For efficiency, we should prerequest all
218  // the data we will need to build the
219  // linear system before doing an element loop.
220  FEBase* elem_fe = NULL;
221  context.get_element_fe(0, elem_fe);
222  elem_fe->get_JxW();
223  elem_fe->get_phi();
224  elem_fe->get_dphi();
225 
226  FEBase* side_fe = NULL;
227  context.get_side_fe( 0, side_fe );
228  side_fe->get_JxW();
229  side_fe->get_phi();
230 
231  FEBase* neighbor_side_fe = NULL;
232  context.get_neighbor_side_fe(0, neighbor_side_fe);
233  neighbor_side_fe->get_phi();
234  }
235 
236  for (const auto & elem : mesh.active_local_element_ptr_range())
237  {
238  context.pre_fe_reinit(system, elem);
239  context.elem_fe_reinit();
240 
241  // Element interior assembly
242  {
243  FEBase* elem_fe = NULL;
244  context.get_element_fe(0, elem_fe);
245 
246  const std::vector<Real> &JxW = elem_fe->get_JxW();
247  const std::vector<std::vector<Real> >& phi = elem_fe->get_phi();
248  const std::vector<std::vector<RealGradient> >& dphi = elem_fe->get_dphi();
249 
250  unsigned int n_dofs = context.get_dof_indices(0).size();
251  unsigned int n_qpoints = context.get_element_qrule().n_points();
252 
253  for (unsigned int qp=0; qp != n_qpoints; qp++)
254  for (unsigned int i=0; i != n_dofs; i++)
255  for (unsigned int j=0; j != n_dofs; j++)
256  context.get_elem_jacobian()(i,j) += JxW[qp] * dphi[i][qp]*dphi[j][qp];
257 
258  for (unsigned int qp=0; qp != n_qpoints; qp++)
259  for (unsigned int i=0; i != n_dofs; i++)
260  context.get_elem_residual()(i) += JxW[qp] * phi[i][qp];
261  }
262 
263  matrix.add_matrix (context.get_elem_jacobian(), context.get_dof_indices());
264  system.rhs->add_vector (context.get_elem_residual(), context.get_dof_indices());
265 
266  // Element side assembly
267  for (context.side = 0; context.side != elem->n_sides(); ++context.side)
268  {
269  // If there is a neighbor, then we proceed with assembly
270  // that involves elem and neighbor
271  const Elem* neighbor = elem->neighbor_ptr(context.get_side());
272  if(neighbor)
273  {
274  context.side_fe_reinit();
275  context.set_neighbor(*neighbor);
276 
277  // This call initializes neighbor data, and also sets
278  // context.dg_terms_are_active() to true
279  context.neighbor_side_fe_reinit();
280 
281  FEBase* side_fe = NULL;
282  context.get_side_fe(0, side_fe);
283 
284  const std::vector<Real> &JxW_face = side_fe->get_JxW();
285  const std::vector<std::vector<Real> >& phi_face = side_fe->get_phi();
286 
287  FEBase* neighbor_side_fe = NULL;
288  context.get_neighbor_side_fe(0, neighbor_side_fe);
289 
290  // These shape functions have been evaluated on the quadrature points
291  // for elem->side on the neighbor element
292  const std::vector<std::vector<Real> >& phi_neighbor_face =
293  neighbor_side_fe->get_phi();
294 
295  const unsigned int n_dofs = context.get_dof_indices(0).size();
296  const unsigned int n_neighbor_dofs = context.get_neighbor_dof_indices(0).size();
297  unsigned int n_sidepoints = context.get_side_qrule().n_points();
298 
299  for (unsigned int qp=0; qp<n_sidepoints; qp++)
300  {
301  for (unsigned int i=0; i<n_dofs; i++)
302  for (unsigned int j=0; j<n_dofs; j++)
303  {
304  context.get_elem_elem_jacobian()(i,j) +=
305  JxW_face[qp] * phi_face[i][qp] * phi_face[j][qp];
306  }
307 
308  for (unsigned int i=0; i<n_dofs; i++)
309  for (unsigned int j=0; j<n_neighbor_dofs; j++)
310  {
311  context.get_elem_neighbor_jacobian()(i,j) +=
312  JxW_face[qp] * phi_face[i][qp] * phi_neighbor_face[j][qp];
313  }
314 
315  for (unsigned int i=0; i<n_neighbor_dofs; i++)
316  for (unsigned int j=0; j<n_neighbor_dofs; j++)
317  {
318  context.get_neighbor_neighbor_jacobian()(i,j) +=
319  JxW_face[qp] * phi_neighbor_face[i][qp] * phi_neighbor_face[j][qp];
320  }
321 
322  for (unsigned int i=0; i<n_neighbor_dofs; i++)
323  for (unsigned int j=0; j<n_dofs; j++)
324  {
325  context.get_neighbor_elem_jacobian()(i,j) +=
326  JxW_face[qp] * phi_neighbor_face[i][qp] * phi_face[j][qp];
327  }
328  }
329 
330  matrix.add_matrix (context.get_elem_elem_jacobian(),
331  context.get_dof_indices(),
332  context.get_dof_indices());
333 
334  matrix.add_matrix (context.get_elem_neighbor_jacobian(),
335  context.get_dof_indices(),
336  context.get_neighbor_dof_indices());
337 
338  matrix.add_matrix (context.get_neighbor_elem_jacobian(),
339  context.get_neighbor_dof_indices(),
340  context.get_dof_indices());
341 
342  matrix.add_matrix (context.get_neighbor_neighbor_jacobian(),
343  context.get_neighbor_dof_indices(),
344  context.get_neighbor_dof_indices());
345  }
346  }
347  }
348 }
Manages consistently variables, degrees of freedom, coefficient vectors, matrices and linear solvers ...
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]...
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
MeshBase & mesh
NumericVector< Number > * rhs
The system matrix.
const T_sys & get_system(std::string_view name) const
This is the MeshBase class.
Definition: mesh_base.h:75
const std::vector< std::vector< OutputGradient > > & get_dphi() const
Definition: fe_base.h:230
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.
This class extends FEMContext in order to provide extra data required to perform local element residu...
virtual_for_inffe const std::vector< Real > & get_JxW() const
Definition: fe_abstract.h:303
void set_neighbor(const unsigned int i, Elem *n)
Assigns n as the neighbor.
Definition: elem.h:2618
const Elem * neighbor_ptr(unsigned int i) const
Definition: elem.h:2598
const MeshBase & get_mesh() const
const SparseMatrix< Number > & get_system_matrix() const
This class forms the foundation from which generic finite elements may be derived.
const std::vector< std::vector< OutputShape > > & get_phi() const
Definition: fe_base.h:207

◆ CPPUNIT_TEST_SUITE_REGISTRATION()

CPPUNIT_TEST_SUITE_REGISTRATION ( SystemsTest  )

◆ cubic_test()

Number cubic_test ( const Point p,
const Parameters ,
const std::string &  ,
const std::string &   
)

Definition at line 351 of file systems_test.C.

References libMesh::Real.

Referenced by TripleFunction::component(), TripleFunction::operator()(), SystemsTest::testProjectCubeWithMeshFunction(), and SystemsTest::tripleValueTest().

355 {
356  const Real & x = p(0);
357  const Real & y = LIBMESH_DIM > 1 ? p(1) : 0;
358  const Real & z = LIBMESH_DIM > 2 ? p(2) : 0;
359 
360  return x*(1-x)*(1-x) + x*x*(1-y) + x*(1-y)*(1-z) + y*(1-y)*z + z*(1-z)*(1-z);
361 }
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real

◆ disc_thirds_test()

Number disc_thirds_test ( const Point p,
const Parameters ,
const std::string &  ,
const std::string &   
)

Definition at line 377 of file systems_test.C.

References libMesh::Real.

Referenced by TripleFunction::component(), TripleFunction::operator()(), and SystemsTest::tripleValueTest().

381 {
382  const Real & x = p(0);
383  const Real & y = LIBMESH_DIM > 1 ? p(1) : 0;
384  const Real & z = LIBMESH_DIM > 2 ? p(2) : 0;
385 
386  return (3*x < 1) + (3*y < 2) + (3*z > 2);
387 }
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real

◆ new_linear_test()

Number new_linear_test ( const Point p,
const Parameters ,
const std::string &  ,
const std::string &   
)

Definition at line 364 of file systems_test.C.

References libMesh::Real.

Referenced by TripleFunction::component(), TripleFunction::operator()(), and SystemsTest::tripleValueTest().

368 {
369  const Real & x = p(0);
370  const Real & y = LIBMESH_DIM > 1 ? p(1) : 0;
371  const Real & z = LIBMESH_DIM > 2 ? p(2) : 0;
372 
373  return x + 2*y + 3*z - 1;
374 }
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real