RE: [Libmesh-users] Periodic boundary conditions

SourceForge Headquarters 225 Broadway Suite 1600 San Diego, CA 92101 +1 (858) 454-5900

The easiest way to find a DOF index from a node number is to use the node
itself.  The Nodes in libMesh are derived from a DofObject, which handles
all of the indexing.

For example:

// get a reference, this is cheaper
// than copying the node
const Node& node = mesh.node(n);

node.dof_number(sys,var,comp);

returns the global degree-of-freedom index for the comp-component (==0 for
lagrange elements) of the var-variable in the sys-system.

Seems to me that this warrants an example in the code so the next person who
wants to solve this problem doesn't have to re-invent the wheel?

-Ben

-----Original Message-----
From: lib...@li...
[mailto:lib...@li...] On Behalf Of Michael
Povolotskyi
Sent: Friday, March 18, 2005 12:31 PM
To: lib...@li...
Subject: [Libmesh-users] Periodic boundary conditions

Dear Ben and Roy,
with your very valuable help I menaged to apply periodic boundary 
conditions with the use
of the DOF constraints.

I post here the function I programmed. It may look not very rational, 
because I, probably, didn't use
Libmesh tools correctly. For example, I didn't find a direct way to find 
a DOF index by a node number,
so I had to loop over all the elements. I should be very happy to 
receive your comments.

Michael.

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

  // It is a good idea to make sure we are assembling
  // the proper system.
  assert (system_name == "Poisson");

  // Declare a performance log.  Give it a descriptive
  // string to identify what part of the code we are
  // logging, since there may be many PerfLogs in an
  // application.
  PerfLog perf_log ("Periodic bc. Assembly",false);

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

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

  // Get a reference to the LinearImplicitSystem we are solving
  LinearImplicitSystem& system = es.get_system<LinearImplicitSystem>
("Poisson");

  // 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();
  DofMap& dof_map = system.get_dof_map();
  // Get a constant reference to the Finite Element type
  // for the first (and only) variable in the system.
  FEType fe_type = dof_map.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 an AutoPtr<FEBase>.  This can be thought
  // of as a pointer that will clean up after itself.  Example 4
  // describes some advantages of  AutoPtr's in the context of
  // quadrature rules.
  AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));

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

  // 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<unsigned int> dof_indices;

  const double pos_tol = 1e-10;
  const double func_tol = 1e-10;

  for (int i = 0; i < mesh.mesh_dimension(); i++) //Loop over all the mesh
directions
    { 
      if  (periodicity[i]) //Check if the periodic b.c. are applied along
the direction i

	{
	  for (unsigned int n = 0; n < mesh.n_nodes(); n++) // Loop over all
the nodes
	    {
	      const Node node1=mesh.node(n);
	      if ( std::abs(node1(i) - min_coord[i]) < pos_tol) // if a node
lies at the boundary
		{ 

		  //let us find dof for it-----------------
		  //
		  unsigned int :: n_dof;
		  bool :: not_found_yet = true;

		  MeshBase::const_element_iterator       el2     =
mesh.active_elements_begin();
		  const MeshBase::const_element_iterator end_el2 =
mesh.active_elements_end();

		  for ( ; ((el2 != end_el2) && not_found_yet)  ; ++el2)  
		    {
		      const Elem* elem = *el2;
		      for (unsigned int i1 = 0; ((i1 < elem->n_nodes()) &&
not_found_yet); i1++)
			{
			  if (node1 == *(elem->get_node(i1)) )
			    {
			      dof_map.dof_indices (elem, dof_indices);
			      n_dof = dof_indices[i1];
			      not_found_yet = false;
			    }
			}
		    }
		  assert(!not_found_yet);

		  //dof is found-------------------------------

		  //let us make a  point that lies at the opposite side

		  Point point2(node1);

		  point2(i) = max_coord[i];

		  //corresponding point is made
                  //--------------------------------------------

		  //let us find an element this point belongs to and
calculate the constraints

		  MeshBase::const_element_iterator       el1     =
mesh.active_elements_begin();
		  const MeshBase::const_element_iterator end_el1 =
mesh.active_elements_end();
		  for ( ; el1 != end_el1 ; ++el1)  
		    {
		      const Elem* elem = *el1;
		      if (elem->contains_point(point2))
			{

			  //elem contains the opposite  point
			  DofMap::DofConstraintRow constraint;

			  dof_map.dof_indices (elem, dof_indices);

			  std::vector<Point> point2_vec(1);

			  point2_vec[0] = point2;

			  std::vector<Point> point2_ref_vec(1);

			  FEInterface::inverse_map (elem->dim(), fe_type ,
elem,  point2_vec,  point2_ref_vec)  ;

			  fe->reinit (elem, &point2_ref_vec);

			  for (int i1 = 0; i1 < phi.size(); i1++)
			    {
			      const std::vector<Point >& qface_point =
fe->get_xyz();
			      if ( std::abs(phi[i1][0]) >  func_tol )
constraint[dof_indices[i1]] = phi[i1][0];
			    }

			  dof_map.add_constraint_row (n_dof,  constraint);

			  break;
			}
		    }

		}

	    }
	}
    }  

}

-------------------------------------------------------
SF email is sponsored by - The IT Product Guide
Read honest & candid reviews on hundreds of IT Products from real users.
Discover which products truly live up to the hype. Start reading now.
http://ads.osdn.com/?ad_id=6595&alloc_id=14396&op=click
_______________________________________________
Libmesh-users mailing list
Lib...@li...
https://lists.sourceforge.net/lists/listinfo/libmesh-users