From: Michael P. <pov...@in...> - 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. |
From: Michael P. <pov...@in...> - 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<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; } } } } } } } |
From: Roy S. <roy...@ic...> - 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 |
From: Michael P. <pov...@in...> - 2005-03-18 19:19:06
|
<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"> <html> <head> <meta http-equiv="Content-Type" content="text/html;charset=US-ASCII"> <title></title> </head> <body text="#000000" bgcolor="#cccccc"> Thank you very much for the hints.<br> I'll try to implement them.<br> Yes, I tested the code and seems it works.<br> Michael.<br> <br> <br> Roy Stogner wrote:<br> <blockquote type="cite" cite="midPine.LNX.4.62.0503181255390.30271@mycroft.localdomain">On Fri, 18 Mar 2005, Michael Povolotskyi wrote: <br> <br> <blockquote type="cite">For example, I didn't find a direct way to find a DOF index by a <br> node number, so I had to loop over all the elements. <br> </blockquote> <br> In general there is no way to find a single DoF index by a node <br> number, since only in the quadratic Lagrange element case do you get a <br> 1-to-1 correspondence between nodes and DoFs. For linear elements on <br> some meshes you get nodes with 0 DoFs, for Hierarchic, Clough-Tocher, <br> etc. elements you get some nodes with multiple DoFs. For adaptive <br> non-Lagrange elements, hanging nodes can have different (but mutually <br> constrained) DoFs depending on which element you're looking at them <br> from! <br> <br> <blockquote type="cite"> for (unsigned int n = 0; n < mesh.n_nodes(); n++) // Loop over all the nodes <br> </blockquote> <br> Rather than looping over every node in the mesh, then every element in <br> the mesh to find out where it is (O(n^2) time cost), you want to try <br> looping over every element in the mesh, then every boundary node on <br> an element (O(n*p^d) time cost, so usually much cheaper). <br> <br> <blockquote type="cite"><br> Point point2(node1); <br> point2(i) = max_coord[i]; <br> //corresponding point is made <br> //-------------------------------------------- <br> //let us find an element this point belongs to and calculate the constraints <br> <br> MeshBase::const_element_iterator el1 = mesh.active_elements_begin(); <br> const MeshBase::const_element_iterator end_el1 = mesh.active_elements_end(); <br> </blockquote> <br> You ought to be able to get more efficient than this by looping over <br> just the coarsest elements, then looping over just the children of <br> whichever coarse element contains your point and recursing until you <br> hit an active element. That would be O(n log n) instead of O(n^2), <br> assuming you start with a very coarse mesh and refine a lot (rather <br> than just starting with a very fine mesh). You could also get to O(n) <br> by making a lookup table (pairing elements with all their periodic <br> neighbors) just once and using it inside your loop. <br> <br> Efficiency aside, this looks like it might work. Have you tested it <br> yet? <br> --- <br> Roy Stogner <br> <br> </blockquote> <br> <br> <pre class="moz-signature" cols="250"> </pre> </body> </html> |
From: Michael P. <pov...@in...> - 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. |
From: Roy S. <roy...@ic...> - 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 |
From: John P. <pet...@cf...> - 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 |
From: Roy S. <roy...@ic...> - 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 |
From: Martin L. <lu...@va...> - 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 |