Work at SourceForge, help us to make it a better place! We have an immediate need for a Support Technician in our San Francisco or Denver office.

## [Libmesh-users] Periodic boundary conditions

 [Libmesh-users] Periodic boundary conditions From: Martin Luethi - 2012-01-02 17:09:48 ```Hi all How are periodic boundary conditions supposed to be used? Example 24 (from Libmesh-0.7.2) does *not* show its usage, since there's never a variable added to the periodic bc, and time stepping stops before the warm patch reaches the boundary. If we add a few time steps, there are ripples in the solution, before it vanishes completely. Ah, OK, I see that the exact_solution is prescribed as a boundary condition, so I think that this example is not helpful at all. Is there a better example for the usage of periodic bcs? Thanks, Martin ```

 [Libmesh-users] Periodic boundary conditions From: Michael Povolotskyi - 2005-03-15 12:53:11 ```Hello everybody, I'm going to solve Poisson equation in a rectangular domain. I'd like to impose periodic boundary conditions at two opposite sides. In order to do this, I have to impose DOF constrains like u_m = u_n, where m and n are DOF's numbers of the points that lie at the opposite sides. Could you, please, tell me how is it possible to impose such constraints in libmesh. Thank you, Michael. ```
 [Libmesh-users] Periodic boundary conditions From: Michael Povolotskyi - 2005-03-18 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 ("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. 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 fe (FEBase::build(dim, fe_type)); // The element shape functions evaluated at the quadrature points. const std::vector >& 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 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 point2_vec(1); point2_vec[0] = point2; std::vector 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& 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; } } } } } } } ```
 Re: [Libmesh-users] Periodic boundary conditions From: Roy Stogner - 2005-03-18 19:09:06 ```On Fri, 18 Mar 2005, Michael Povolotskyi wrote: > 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. In general there is no way to find a single DoF index by a node number, since only in the quadratic Lagrange element case do you get a 1-to-1 correspondence between nodes and DoFs. For linear elements on some meshes you get nodes with 0 DoFs, for Hierarchic, Clough-Tocher, etc. elements you get some nodes with multiple DoFs. For adaptive non-Lagrange elements, hanging nodes can have different (but mutually constrained) DoFs depending on which element you're looking at them from! > for (unsigned int n = 0; n < mesh.n_nodes(); n++) // Loop over all > the nodes Rather than looping over every node in the mesh, then every element in the mesh to find out where it is (O(n^2) time cost), you want to try looping over every element in the mesh, then every boundary node on an element (O(n*p^d) time cost, so usually much cheaper). > > 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(); You ought to be able to get more efficient than this by looping over just the coarsest elements, then looping over just the children of whichever coarse element contains your point and recursing until you hit an active element. That would be O(n log n) instead of O(n^2), assuming you start with a very coarse mesh and refine a lot (rather than just starting with a very fine mesh). You could also get to O(n) by making a lookup table (pairing elements with all their periodic neighbors) just once and using it inside your loop. Efficiency aside, this looks like it might work. Have you tested it yet? --- Roy Stogner ```
 Re: [Libmesh-users] Periodic boundary conditions From: Michael Povolotskyi - 2005-03-18 19:19:06 ``` Thank you very much for the hints.
I'll try to implement them.
Yes, I tested the code and seems it works.
Michael.

Roy Stogner wrote:
On Fri, 18 Mar 2005, Michael Povolotskyi wrote:

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.

In general there is no way to find a single DoF index by a node
number, since only in the quadratic Lagrange element case do you get a
1-to-1 correspondence between nodes and DoFs.  For linear elements on
some meshes you get nodes with 0 DoFs, for Hierarchic, Clough-Tocher,
etc. elements you get some nodes with multiple DoFs.  For adaptive
non-Lagrange elements, hanging nodes can have different (but mutually
constrained) DoFs depending on which element you're looking at them
from!

for (unsigned int n = 0; n < mesh.n_nodes(); n++) // Loop over all the nodes

Rather than looping over every node in the mesh, then every element in
the mesh to find out where it is (O(n^2) time cost), you want to try
looping over every element in the mesh, then every boundary node on
an element (O(n*p^d) time cost, so usually much cheaper).

Point point2(node1);
point2(i) = max_coord[i];
//--------------------------------------------
//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();

You ought to be able to get more efficient than this by looping over
just the coarsest elements, then looping over just the children of
whichever coarse element contains your point and recursing until you
hit an active element.  That would be O(n log n) instead of O(n^2),
assuming you start with a very coarse mesh and refine a lot (rather
than just starting with a very fine mesh).  You could also get to O(n)
by making a lookup table (pairing elements with all their periodic
neighbors) just once and using it inside your loop.

Efficiency aside, this looks like it might work.  Have you tested it
yet?
---
Roy Stogner

```
 Re: [Libmesh-users] Periodic boundary conditions From: Michael Povolotskyi - 2005-03-21 11:09:56 ```> > You ought to be able to get more efficient than this by looping over > just the coarsest elements, then looping over just the children of > whichever coarse element contains your point and recursing until you > hit an active element. That would be O(n log n) instead of O(n^2), > assuming you start with a very coarse mesh and refine a lot (rather > than just starting with a very fine mesh). Thank you very much for the hint! Does the special iterator exist that allows to loop over the coarsest elements? Thank you, Michael. ```
 Re: [Libmesh-users] Periodic boundary conditions From: Roy Stogner - 2005-03-21 16:05:28 ```On Mon, 21 Mar 2005, Michael Povolotskyi wrote: >> You ought to be able to get more efficient than this by looping over >> just the coarsest elements, then looping over just the children of >> whichever coarse element contains your point and recursing until you >> hit an active element. That would be O(n log n) instead of O(n^2), >> assuming you start with a very coarse mesh and refine a lot (rather >> than just starting with a very fine mesh). > > Thank you very much for the hint! > Does the special iterator exist that allows to loop over the coarsest > elements? Hmmm... the iterator class exists, but the function to create it doesn't seem to. The level_elem_iterator class for the coarsest ought to be created by mesh->level_elements_begin(0), but the only related function I can find is MeshBase::not_level_elements_begin(). John? --- Roy Stogner ```
 Re: [Libmesh-users] Periodic boundary conditions From: John Peterson - 2005-03-21 15:17:21 ```Roy Stogner writes: > On Mon, 21 Mar 2005, Michael Povolotskyi wrote: > > >> You ought to be able to get more efficient than this by looping over > >> just the coarsest elements, then looping over just the children of > >> whichever coarse element contains your point and recursing until you > >> hit an active element. That would be O(n log n) instead of O(n^2), > >> assuming you start with a very coarse mesh and refine a lot (rather > >> than just starting with a very fine mesh). > > > > Thank you very much for the hint! > > Does the special iterator exist that allows to loop over the coarsest > > elements? > > Hmmm... the iterator class exists, but the function to create it > doesn't seem to. The level_elem_iterator class for the coarsest ought > to be created by mesh->level_elements_begin(0), but the only related > function I can find is MeshBase::not_level_elements_begin(). > John? It can be added ASAP today. It was just never used before so it never existed... -J ```
 [Libmesh-users] Periodic boundary conditions From: Roy Stogner - 2007-06-01 01:44:54 ```I've just added support for using DofMap constraints to impose periodic boundary conditions. There's been some interest in this on the mailing lists in the past, so if anyone would like to test it out: Download the latest CVS. (Or just "cvs update" an existing tree and possibly "make distclean" to be safe) Run "./configure --enable-everything" (or just "--enable-periodic" along with whatever other options you need), and recompile. Create a mesh with boundary ids corresponding to each side of the periodic boundary. The coarse mesh needs to be conforming across the boundary. Add code like the following (which works with build_square's boundary ids on a square of width W, height H) before equation_systems.init(): PeriodicBoundary p1, p2; p1.myboundary = 0; // Bottom p1.pairedboundary = 2; // Top p1.translation_vector = RealVectorValue(0, H); p2.myboundary = 3; // Left p2.pairedboundary = 1; // Right p2.translation_vector = RealVectorValue(W, 0); mysystem.get_dof_map().add_periodic_boundary(p1); mysystem.get_dof_map().add_periodic_boundary(p2); In theory, this feature should work even with non-straight boundaries, periodic boundaries that meet at corners, adaptive refinement (except that level-1 rules aren't yet enforced between boundaries), C0 and C1 elements, etc. In practice, I've only tested it on one application on a uniform square grid, and it'll probably break if you look at it funny. Send any bug reports to libmesh-devel. API suggestions are appreciated too. I've currently got the PeriodicBoundaries defined in the DofMap, but one could argue that it should be a property of a Mesh or that it should be per-variable. --- Roy ```
 [Libmesh-users] Periodic boundary conditions From: Martin Luethi - 2012-01-02 17:09:48 ```Hi all How are periodic boundary conditions supposed to be used? Example 24 (from Libmesh-0.7.2) does *not* show its usage, since there's never a variable added to the periodic bc, and time stepping stops before the warm patch reaches the boundary. If we add a few time steps, there are ripples in the solution, before it vanishes completely. Ah, OK, I see that the exact_solution is prescribed as a boundary condition, so I think that this example is not helpful at all. Is there a better example for the usage of periodic bcs? Thanks, Martin ```