amr.C File Reference

Go to the source code of this file.

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

Function Documentation

assemble() [1/2]

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

Definition at line 239 of file miscellaneous_ex4.C.

References libMesh::SparseMatrix< T >::add_matrix(), libMesh::NumericVector< T >::add_vector(), libMesh::FEGenericBase< OutputType >::build(), libMesh::NumericVector< T >::close(), libMesh::SparseMatrix< T >::close(), libMesh::DofMap::constrain_element_dyad_matrix(), libMesh::DofMap::constrain_element_matrix_and_vector(), dim, libMesh::DofMap::dof_indices(), libMesh::System::get_dof_map(), libMesh::System::get_matrix(), libMesh::EquationSystems::get_mesh(), libMesh::EquationSystems::get_system(), libMesh::ImplicitSystem::get_system_matrix(), libMesh::System::get_vector(), libMesh::libmesh_ignore(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::Real, libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), libMesh::ExplicitSystem::rhs, and libMesh::System::variable_type().

Referenced by main().

{
  // Ignore unused parameter warnings when !LIBMESH_ENABLE_AMR.
  libmesh_ignore(es, system_name);

#ifdef LIBMESH_ENABLE_AMR
  // It is a good idea to make sure we are assembling
  // the proper system.
  libmesh_assert_equal_to (system_name, "System");

  // Get a constant reference to the mesh object.
  const MeshBase & mesh = es.get_mesh();

  // The dimension that we are running
  const unsigned int dim = mesh.mesh_dimension();

  // Get a reference to the Convection-Diffusion system object.
  LinearImplicitSystem & system =
    es.get_system<LinearImplicitSystem> ("System");

  // Get the Finite Element type for the first (and only)
  // variable in the system.
  FEType fe_type = system.variable_type(0);

  // Build a Finite Element object of the specified type.  Since the
  // FEBase::build() member dynamically creates memory we will
  // store the object as a std::unique_ptr<FEBase>.  This can be thought
  // of as a pointer that will clean up after itself.
  std::unique_ptr<FEBase> fe      (FEBase::build(dim, fe_type));
  std::unique_ptr<FEBase> fe_face (FEBase::build(dim, fe_type));

  // A Gauss quadrature rule for numerical integration.
  // Let the FEType object decide what order rule is appropriate.
  QGauss qrule (dim,   fe_type.default_quadrature_order());
  QGauss qface (dim-1, fe_type.default_quadrature_order());

  // Tell the finite element object to use our quadrature rule.
  fe->attach_quadrature_rule      (&qrule);
  fe_face->attach_quadrature_rule (&qface);

  // Here we define some references to cell-specific data that
  // will be used to assemble the linear system.  We will start
  // with the element Jacobian * quadrature weight at each integration point.
  const std::vector<Real> & JxW      = fe->get_JxW();
  const std::vector<Real> & JxW_face = fe_face->get_JxW();

  // The element shape functions evaluated at the quadrature points.
  const std::vector<std::vector<Real>> & phi = fe->get_phi();
  const std::vector<std::vector<Real>> & psi = fe_face->get_phi();

  // The element shape function gradients evaluated at the quadrature
  // points.
  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();

  // The XY locations of the quadrature points used for face integration
  //const std::vector<Point>& qface_points = fe_face->get_xyz();

  // A reference to the DofMap object for this system.  The DofMap
  // object handles the index translation from node and element numbers
  // to degree of freedom numbers.  We will talk more about the DofMap
  // in future examples.
  const DofMap & dof_map = system.get_dof_map();

  // Define data structures to contain the element matrix
  // and right-hand-side vector contribution.  Following
  // basic finite element terminology we will denote these
  // "Ke" and "Fe".
  DenseMatrix<Number> Ke;
  DenseVector<Number> Fe;

  // Analogous data structures for thw two vectors v and w that form
  // the tensor shell matrix.
  DenseVector<Number> Ve;
  DenseVector<Number> We;

  // This vector will hold the degree of freedom indices for
  // the element.  These define where in the global system
  // the element degrees of freedom get mapped.
  std::vector<dof_id_type> dof_indices;

  // The global system matrix
  SparseMatrix<Number> & matrix = system.get_system_matrix();

  // Now we will loop over all the elements in the mesh that
  // live on the local processor. We will compute the element
  // matrix and right-hand-side contribution.  Since the mesh
  // will be refined we want to only consider the ACTIVE elements,
  // hence we use a variant of the active_elem_iterator.
  for (const auto & elem : mesh.active_local_element_ptr_range())
    {
      // Get the degree of freedom indices for the
      // current element.  These define where in the global
      // matrix and right-hand-side this element will
      // contribute to.
      dof_map.dof_indices (elem, dof_indices);

      // Compute the element-specific data for the current
      // element.  This involves computing the location of the
      // quadrature points (q_point) and the shape functions
      // (phi, dphi) for the current element.
      fe->reinit (elem);

      // Zero the element matrix and right-hand side before
      // summing them.  We use the resize member here because
      // the number of degrees of freedom might have changed from
      // the last element.  Note that this will be the case if the
      // element type is different (i.e. the last element was a
      // triangle, now we are on a quadrilateral).
      const unsigned int n_dofs =
        cast_int<unsigned int>(dof_indices.size());

      Ke.resize (n_dofs, n_dofs);

      Fe.resize (n_dofs);
      Ve.resize (n_dofs);
      We.resize (n_dofs);

      // Now we will build the element matrix and right-hand-side.
      // Constructing the RHS requires the solution and its
      // gradient from the previous timestep.  This myst be
      // calculated at each quadrature point by summing the
      // solution degree-of-freedom values by the appropriate
      // weight functions.
      for (unsigned int qp=0; qp<qrule.n_points(); qp++)
        {
          // Now compute the element matrix and RHS contributions.
          for (unsigned int i=0; i<n_dofs; i++)
            {
              // The RHS contribution
              Fe(i) += JxW[qp]*phi[i][qp];

              for (unsigned int j=0; j<n_dofs; j++)
                {
                  // The matrix contribution
                  Ke(i,j) += JxW[qp]*(
                                      // Stiffness matrix
                                      (dphi[i][qp]*dphi[j][qp])
                                      );
                }

              // V and W are the same for this example.
              Ve(i) += JxW[qp]*phi[i][qp];
              We(i) += JxW[qp]*phi[i][qp];
            }
        }

      // At this point the interior element integration has
      // been completed.  However, we have not yet addressed
      // boundary conditions.  For this example we will only
      // consider simple Dirichlet boundary conditions imposed
      // via the penalty method.
      //
      // The following loops over the sides of the element.
      // If the element has no neighbor on a side then that
      // side MUST live on a boundary of the domain.
      {
        // The penalty value.
        const Real penalty = 1.e10;

        // The following loops over the sides of the element.
        // If the element has no neighbor on a side then that
        // side MUST live on a boundary of the domain.
        for (auto s : elem->side_index_range())
          if (elem->neighbor_ptr(s) == nullptr)
            {
              fe_face->reinit(elem, s);

              for (unsigned int qp=0; qp<qface.n_points(); qp++)
                {
                  // Matrix contribution
                  for (unsigned int i=0; i<n_dofs; i++)
                    for (unsigned int j=0; j<n_dofs; j++)
                      Ke(i,j) += penalty*JxW_face[qp]*psi[i][qp]*psi[j][qp];
                }
            }
      }


      // We have now built the element matrix and RHS vector in terms
      // of the element degrees of freedom.  However, it is possible
      // that some of the element DOFs are constrained to enforce
      // solution continuity, i.e. they are not really "free".  We need
      // to constrain those DOFs in terms of non-constrained DOFs to
      // ensure a continuous solution.  The
      // DofMap::constrain_element_matrix_and_vector() method does
      // just that.

      // However, constraining both the sparse matrix (and right hand
      // side) plus the rank 1 matrix is tricky.  The dof_indices
      // vector has to be backed up for that because the constraining
      // functions modify it.

      std::vector<dof_id_type> dof_indices_backup(dof_indices);
      dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
      dof_indices = dof_indices_backup;
      dof_map.constrain_element_dyad_matrix(Ve, We, dof_indices);

      // The element matrix and right-hand-side are now built
      // for this element.  Add them to the global matrix and
      // right-hand-side vector.  The SparseMatrix::add_matrix()
      // and NumericVector::add_vector() members do this for us.
      matrix.add_matrix (Ke, dof_indices);
      system.get_matrix("Preconditioner").add_matrix (Ke, dof_indices);
      system.rhs->add_vector (Fe, dof_indices);
      system.get_vector("v").add_vector(Ve, dof_indices);
      system.get_vector("w").add_vector(We, dof_indices);
    }
  // Finished computing the system matrix and right-hand side.

  // Matrices and vectors must be closed manually.  This is necessary
  // because the matrix is not directly used as the system matrix (in
  // which case the solver closes it) but as a part of a shell matrix.
  matrix.close();
  system.get_matrix("Preconditioner").close();
  system.rhs->close();
  system.get_vector("v").close();
  system.get_vector("w").close();

#endif // #ifdef LIBMESH_ENABLE_AMR
}

assemble() [2/2]

void assemble	(	EquationSystems &	es,
		const std::string &	libmesh_dbg_varsystem_name
	)

Definition at line 135 of file amr.C.

References libMesh::SparseMatrix< T >::add_matrix(), libMesh::NumericVector< T >::add_vector(), libMesh::FEGenericBase< OutputType >::build(), libMesh::DofMap::constrain_element_matrix_and_vector(), dim, libMesh::DofMap::dof_indices(), libMesh::FIFTH, libMesh::FIRST, libMesh::System::get_dof_map(), libMesh::EquationSystems::get_mesh(), libMesh::EquationSystems::get_system(), libMesh::ImplicitSystem::get_system_matrix(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::QBase::n_points(), libMesh::out, libMesh::Real, libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), and libMesh::ExplicitSystem::rhs.

{
  libmesh_assert_equal_to (system_name, "primary");

  const MeshBase & mesh   = es.get_mesh();
  const unsigned int dim = mesh.mesh_dimension();

  // Also use a 3x3x3 quadrature rule (3D).  Then tell the FE
  // about the geometry of the problem and the quadrature rule
  FEType fe_type (FIRST);

  std::unique_ptr<FEBase> fe(FEBase::build(dim, fe_type));
  QGauss qrule(dim, FIFTH);

  fe->attach_quadrature_rule (&qrule);

  std::unique_ptr<FEBase> fe_face(FEBase::build(dim, fe_type));
  QGauss qface(dim-1, FIFTH);

  fe_face->attach_quadrature_rule(&qface);

  LinearImplicitSystem & system =
    es.get_system<LinearImplicitSystem>("primary");


  // These are references to cell-specific data
  const std::vector<Real> & JxW_face                   = fe_face->get_JxW();
  const std::vector<Real> & JxW                        = fe->get_JxW();
  const std::vector<Point> & q_point                   = fe->get_xyz();
  const std::vector<std::vector<Real>> & phi          = fe->get_phi();
  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();

  std::vector<dof_id_type> dof_indices_U;
  std::vector<dof_id_type> dof_indices_V;
  const DofMap & dof_map = system.get_dof_map();

  DenseMatrix<Number> Kuu;
  DenseMatrix<Number> Kvv;
  DenseVector<Number> Fu;
  DenseVector<Number> Fv;

  Real vol=0., area=0.;

  SparseMatrix<Number> & matrix = system.get_system_matrix();

  for (const auto & elem : mesh.active_local_element_ptr_range())
    {
      // recompute the element-specific data for the current element
      fe->reinit (elem);


      //fe->print_info();

      dof_map.dof_indices(elem, dof_indices_U, 0);
      dof_map.dof_indices(elem, dof_indices_V, 1);

      const unsigned int n_phi = cast_int<unsigned int>(phi.size());

      // zero the element matrix and vector
      Kuu.resize (n_phi, n_phi);

      Kvv.resize (n_phi, n_phi);

      Fu.resize (n_phi);
      Fv.resize (n_phi);

      // standard stuff...  like in code 1.
      for (unsigned int gp=0; gp<qrule.n_points(); gp++)
        {
          for (unsigned int i=0; i<n_phi; ++i)
            {
              // this is tricky.  ig is the _global_ dof index corresponding
              // to the _global_ vertex number elem->node_id(i).  Note that
              // in general these numbers will not be the same (except for
              // the case of one unknown per node on one subdomain) so
              // we need to go through the dof_map

              const Real f = q_point[gp]*q_point[gp];
              //    const Real f = (q_point[gp](0) +
              //    q_point[gp](1) +
              //    q_point[gp](2));

              // add jac*weight*f*phi to the RHS in position ig

              Fu(i) += JxW[gp]*f*phi[i][gp];
              Fv(i) += JxW[gp]*f*phi[i][gp];

              for (unsigned int j=0; j != n_phi; ++j)
                {

                  Kuu(i,j) += JxW[gp]*((phi[i][gp])*(phi[j][gp]));

                  Kvv(i,j) += JxW[gp]*((phi[i][gp])*(phi[j][gp]) +
                                       1.*((dphi[i][gp])*(dphi[j][gp])));
                };
            };
          vol += JxW[gp];
        };


      // You can't compute "area" (perimeter) if you are in 2D
      if (dim == 3)
        {
          for (auto side : elem->side_index_range())
            if (elem->neighbor_ptr(side) == nullptr)
              {
                fe_face->reinit (elem, side);

                for (const auto & val : JxW_face)
                  area += val;
              }
        }

      // Constrain the DOF indices.
      dof_map.constrain_element_matrix_and_vector(Kuu, Fu, dof_indices_U);
      dof_map.constrain_element_matrix_and_vector(Kvv, Fv, dof_indices_V);


      system.rhs->add_vector(Fu, dof_indices_U);
      system.rhs->add_vector(Fv, dof_indices_V);

      matrix.add_matrix(Kuu, dof_indices_U);
      matrix.add_matrix(Kvv, dof_indices_V);
    }

  libMesh::out << "Vol="  << vol << std::endl;

  if (dim == 3)
    libMesh::out << "Area=" << area << std::endl;
}

main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 46 of file amr.C.

References libMesh::DofMap::_dof_coupling, libMesh::EquationSystems::add_system(), libMesh::System::add_variable(), assemble(), libMesh::System::attach_assemble_function(), dim, libMesh::MeshBase::elem_ref(), libMesh::FIRST, libMesh::System::get_dof_map(), libMesh::TriangleWrapper::init(), libMesh::EquationSystems::init(), mesh, libMesh::DofMap::print_dof_constraints(), libMesh::EquationSystems::print_info(), libMesh::MeshBase::print_info(), libMesh::MeshBase::read(), libMesh::Elem::REFINE, libMesh::MeshRefinement::refine_and_coarsen_elements(), libMesh::EquationSystems::reinit(), libMesh::CouplingMatrix::resize(), libMesh::Elem::set_refinement_flag(), libMesh::LinearImplicitSystem::solve(), libMesh::MeshRefinement::uniformly_refine(), libMesh::GMVIO::write(), and libMesh::MeshOutput< MT >::write_equation_systems().

{
  LibMeshInit init(argc, argv);

  if (argc < 4)
    {
      std::cout << "Usage: ./prog -d DIM filename" << std::endl;
      libmesh_terminate();
    }

  // Variables to get us started
  const unsigned char dim = cast_int<unsigned char>(atoi(argv[2]));

  std::string meshname  (argv[3]);

  // declare a mesh...
  Mesh mesh(init.comm(), dim);

  // Read a mesh
  mesh.read(meshname);

  GMVIO(mesh).write ("out_0.gmv");

  mesh.elem_ref(0).set_refinement_flag (Elem::REFINE);

  MeshRefinement mesh_refinement (mesh);

  mesh_refinement.refine_and_coarsen_elements ();
  mesh_refinement.uniformly_refine (2);

  mesh.print_info();


  // Set up the equation system(s)
  EquationSystems es (mesh);

  LinearImplicitSystem & primary =
    es.add_system<LinearImplicitSystem>("primary");

  primary.add_variable ("U", FIRST);
  primary.add_variable ("V", FIRST);

  primary.get_dof_map()._dof_coupling->resize(2);
  (*primary.get_dof_map()._dof_coupling)(0,0) = 1;
  (*primary.get_dof_map()._dof_coupling)(1,1) = 1;

  primary.attach_assemble_function(assemble);

  es.init ();

  es.print_info ();
  primary.get_dof_map().print_dof_constraints ();

  // call the solver.
  primary.solve ();

  GMVIO(mesh).write_equation_systems ("out_1.gmv",
                                      es);



  // Refine uniformly
  mesh_refinement.uniformly_refine (1);
  es.reinit ();

  // Write out the projected solution
  GMVIO(mesh).write_equation_systems ("out_2.gmv",
                                      es);

  // Solve again. Output the refined solution
  primary.solve ();
  GMVIO(mesh).write_equation_systems ("out_3.gmv",
                                      es);

  return 0;
}