From: Michael Povolotskyi <povolotskyi@in...>  20050318 18:31:41

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; } } } } } } } 