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 = 1e10;
const double func_tol = 1e10;
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;
}
}
}
}
}
}
}
