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 &system_name)
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 &	system_name
	)

Definition at line 109 of file systems_test.C.

{
  const MeshBase& mesh = es.get_mesh();
  LinearImplicitSystem& system = es.get_system<LinearImplicitSystem>("test");
  const DofMap& dof_map = system.get_dof_map();
 
  DenseMatrix<Number> Ke;
  DenseVector<Number> Fe;
 
  std::vector<dof_id_type> dof_indices;
 
  MeshBase::const_element_iterator       el     = mesh.active_local_elements_begin();
  const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end();
 
  for ( ; el != end_el; ++el)
    {
      const Elem* elem = *el;
 
      if(elem->type() == NODEELEM)
      {
        continue;
      }
 
      dof_map.dof_indices (elem, dof_indices);
      const unsigned int n_dofs = dof_indices.size();
 
      Ke.resize (n_dofs, n_dofs);
      Fe.resize (n_dofs);
 
      for(unsigned int i=0; i<n_dofs; i++)
      {
        Ke(i,i) = 1.;
        Fe(i) = 1.;
      }
 
      system.matrix->add_matrix (Ke, dof_indices);
      system.rhs->add_vector    (Fe, dof_indices);
    }
 
  // Add matrix for extra coupled dofs
  {
    const Node & node_1 = mesh.node_ref(1);
    const Node & node_2 = mesh.node_ref(2);
    dof_indices.resize(6);
    dof_indices[0] =
      node_1.dof_number(system.number(), system.variable_number("u"), 0);
    dof_indices[1] =
      node_1.dof_number(system.number(), system.variable_number("v"), 0);
    dof_indices[2] =
      node_1.dof_number(system.number(), system.variable_number("w"), 0);
 
    dof_indices[3] =
      node_2.dof_number(system.number(), system.variable_number("u"), 0);
    dof_indices[4] =
      node_2.dof_number(system.number(), system.variable_number("v"), 0);
    dof_indices[5] =
      node_2.dof_number(system.number(), system.variable_number("w"), 0);
 
    const unsigned int n_dofs = dof_indices.size();
    Ke.resize (n_dofs, n_dofs);
    Fe.resize (n_dofs);
 
    for(unsigned int i=0; i<n_dofs; i++)
    {
      Ke(i,i) = 1.;
      Fe(i) = 1.;
    }
 
    system.matrix->add_matrix (Ke, dof_indices);
    system.rhs->add_vector    (Fe, dof_indices);
  }
 
  system.rhs->close();
  system.matrix->close();
}

References libMesh::MeshBase::active_local_elements_begin(), libMesh::MeshBase::active_local_elements_end(), libMesh::SparseMatrix< T >::add_matrix(), libMesh::NumericVector< T >::add_vector(), libMesh::SparseMatrix< T >::close(), libMesh::NumericVector< T >::close(), libMesh::DofObject::dof_number(), libMesh::System::get_dof_map(), libMesh::EquationSystems::get_mesh(), libMesh::EquationSystems::get_system(), libMesh::ImplicitSystem::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 SystemsTest::testDofCouplingWithVarGroups().

assembly_with_dg_fem_context()

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

Definition at line 187 of file systems_test.C.

{
  const MeshBase& mesh = es.get_mesh();
  LinearImplicitSystem& system = es.get_system<LinearImplicitSystem>("test");
 
  DenseMatrix<Number> Ke;
  DenseVector<Number> Fe;
 
  std::vector<dof_id_type> dof_indices;
 
  DGFEMContext context(system);
  {
    // For efficiency, we should prerequest all
    // the data we will need to build the
    // linear system before doing an element loop.
    FEBase* elem_fe = NULL;
    context.get_element_fe(0, elem_fe);
    elem_fe->get_JxW();
    elem_fe->get_phi();
    elem_fe->get_dphi();
 
    FEBase* side_fe = NULL;
    context.get_side_fe( 0, side_fe );
    side_fe->get_JxW();
    side_fe->get_phi();
  }
 
  for (const auto & elem : mesh.active_local_element_ptr_range())
    {
      context.pre_fe_reinit(system, elem);
      context.elem_fe_reinit();
 
      // Element interior assembly
      {
        FEBase* elem_fe = NULL;
        context.get_element_fe(0, elem_fe);
 
        const std::vector<Real> &JxW = elem_fe->get_JxW();
        const std::vector<std::vector<Real> >& phi = elem_fe->get_phi();
        const std::vector<std::vector<RealGradient> >& dphi = elem_fe->get_dphi();
 
        unsigned int n_dofs = context.get_dof_indices(0).size();
        unsigned int n_qpoints = context.get_element_qrule().n_points();
 
        for (unsigned int qp=0; qp != n_qpoints; qp++)
          for (unsigned int i=0; i != n_dofs; i++)
            for (unsigned int j=0; j != n_dofs; j++)
              context.get_elem_jacobian()(i,j) += JxW[qp] * dphi[i][qp]*dphi[j][qp];
 
        for (unsigned int qp=0; qp != n_qpoints; qp++)
          for (unsigned int i=0; i != n_dofs; i++)
            context.get_elem_residual()(i) += JxW[qp] * phi[i][qp];
      }
 
      system.matrix->add_matrix (context.get_elem_jacobian(), context.get_dof_indices());
      system.rhs->add_vector (context.get_elem_residual(), context.get_dof_indices());
 
      // Element side assembly
      for (context.side = 0; context.side != elem->n_sides(); ++context.side)
        {
          // If there is a neighbor, then we proceed with assembly
          // that involves elem and neighbor
          const Elem* neighbor = elem->neighbor_ptr(context.get_side());
          if(neighbor)
            {
              context.side_fe_reinit();
              context.set_neighbor(*neighbor);
 
              // This call initializes neighbor data, and also sets
              // context.dg_terms_are_active() to true
              context.neighbor_side_fe_reinit();
 
              FEBase* side_fe = NULL;
              context.get_side_fe(0, side_fe);
 
              const std::vector<Real> &JxW_face = side_fe->get_JxW();
              const std::vector<std::vector<Real> >& phi_face = side_fe->get_phi();
 
              FEBase* neighbor_side_fe = NULL;
              context.get_neighbor_side_fe(0, neighbor_side_fe);
 
              // These shape functions have been evaluated on the quadrature points
              // for elem->side on the neighbor element
              const std::vector<std::vector<Real> >& phi_neighbor_face =
                neighbor_side_fe->get_phi();
 
              const unsigned int n_dofs = context.get_dof_indices(0).size();
              const unsigned int n_neighbor_dofs = context.get_neighbor_dof_indices(0).size();
              unsigned int n_sidepoints = context.get_side_qrule().n_points();
 
              for (unsigned int qp=0; qp<n_sidepoints; qp++)
                {
                  for (unsigned int i=0; i<n_dofs; i++)
                    for (unsigned int j=0; j<n_dofs; j++)
                      {
                        context.get_elem_elem_jacobian()(i,j) +=
                          JxW_face[qp] * phi_face[i][qp] * phi_face[j][qp];
                      }
 
                  for (unsigned int i=0; i<n_dofs; i++)
                    for (unsigned int j=0; j<n_neighbor_dofs; j++)
                      {
                        context.get_elem_neighbor_jacobian()(i,j) +=
                          JxW_face[qp] * phi_face[i][qp] * phi_neighbor_face[j][qp];
                      }
 
                  for (unsigned int i=0; i<n_neighbor_dofs; i++)
                    for (unsigned int j=0; j<n_neighbor_dofs; j++)
                      {
                        context.get_neighbor_neighbor_jacobian()(i,j) +=
                          JxW_face[qp] * phi_neighbor_face[i][qp] * phi_neighbor_face[j][qp];
                      }
 
                  for (unsigned int i=0; i<n_neighbor_dofs; i++)
                    for (unsigned int j=0; j<n_dofs; j++)
                      {
                        context.get_neighbor_elem_jacobian()(i,j) +=
                          JxW_face[qp] * phi_neighbor_face[i][qp] * phi_face[j][qp];
                      }
                }
 
              system.matrix->add_matrix (context.get_elem_elem_jacobian(),
                                         context.get_dof_indices(),
                                         context.get_dof_indices());
 
              system.matrix->add_matrix (context.get_elem_neighbor_jacobian(),
                                         context.get_dof_indices(),
                                         context.get_neighbor_dof_indices());
 
              system.matrix->add_matrix (context.get_neighbor_elem_jacobian(),
                                         context.get_neighbor_dof_indices(),
                                         context.get_dof_indices());
 
              system.matrix->add_matrix (context.get_neighbor_neighbor_jacobian(),
                                         context.get_neighbor_dof_indices(),
                                         context.get_neighbor_dof_indices());
            }
        }
      }
}

References libMesh::MeshBase::active_local_element_ptr_range(), 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::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().

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 330 of file systems_test.C.

{
  const Real & x = p(0);
  const Real & y = LIBMESH_DIM > 1 ? p(1) : 0;
  const Real & z = LIBMESH_DIM > 2 ? p(2) : 0;
 
  return x*(1-x)*(1-x) + x*x*(1-y) + x*(1-y)*(1-z) + y*(1-y)*z + z*(1-z)*(1-z);
}

References libMesh::Real.

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

disc_thirds_test()

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

Definition at line 356 of file systems_test.C.

{
  const Real & x = p(0);
  const Real & y = LIBMESH_DIM > 1 ? p(1) : 0;
  const Real & z = LIBMESH_DIM > 2 ? p(2) : 0;
 
  return (3*x < 1) + (3*y < 2) + (3*z > 2);
}

References libMesh::Real.

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

new_linear_test()

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

Definition at line 343 of file systems_test.C.

{
  const Real & x = p(0);
  const Real & y = LIBMESH_DIM > 1 ? p(1) : 0;
  const Real & z = LIBMESH_DIM > 2 ? p(2) : 0;
 
  return x + 2*y + 3*z - 1;
}

References libMesh::Real.

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