From: KIRK, BENJAMIN (JSCEG) (NASA) <benjamin.kirk1@na...>  20050315 14:19:44

If you know the degree of freedom indices it should be easy. Essentially what you want to do is replace one set of degrees of freedom with another. Assuming u_m corresponds to DOFs n_u_m and similarly u_n is n_u_n you should be able to do something like this: DofConstraintRow constraint; constraint[n_u_n] = 1; dof_map.add_constraint_row (n_u_m, constraint); Essentially this imposes u_m = u_n. Please let me know how this works out. Ben Original Message From: libmeshusersadmin@... [mailto:libmeshusersadmin@...] On Behalf Of Michael Povolotskyi Sent: Tuesday, March 15, 2005 6:52 AM To: libmesh  ML Subject: [Libmeshusers] Periodic boundary conditions 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.  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 _______________________________________________ Libmeshusers mailing list Libmeshusers@... https://lists.sourceforge.net/lists/listinfo/libmeshusers 
From: KIRK, BENJAMIN (JSCEG) (NASA) <benjamin.kirk1@na...>  20050318 19:27:00

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 degreeoffreedom index for the compcomponent (==0 for lagrange elements) of the varvariable in the syssystem. 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 reinvent the wheel? Ben Original Message From: libmeshusersadmin@... [mailto:libmeshusersadmin@...] On Behalf Of Michael Povolotskyi Sent: Friday, March 18, 2005 12:31 PM To: libmeshusers@... Subject: [Libmeshusers] 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 = 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; } } } } } } }  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 _______________________________________________ Libmeshusers mailing list Libmeshusers@... https://lists.sourceforge.net/lists/listinfo/libmeshusers 
From: Roy Stogner <roystgnr@ic...>  20050318 19:36:16

On Fri, 18 Mar 2005, KIRK, BENJAMIN (JSCEG) (NASA) wrote: > 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 reinvent the wheel? I hope not. I'm making some headway getting Lagrangespecific functions extended to handle more general elements; I'd rather not encourage the creation of more Lagrangeonly code!  Roy Stogner 
From: KIRK, BENJAMIN (JSCEG) (NASA) <benjamin.kirk1@na...>  20050318 19:38:22

Of course with proper generalization... :) BTW, are there any difficulties removing those insidious Lagrange assumptions you need help with? Ben Original Message From: Roy Stogner [mailto:roystgnr@...] Sent: Friday, March 18, 2005 1:36 PM To: KIRK, BENJAMIN (JSCEG) (NASA) Cc: libmeshusers@... Subject: RE: [Libmeshusers] Periodic boundary conditions On Fri, 18 Mar 2005, KIRK, BENJAMIN (JSCEG) (NASA) wrote: > 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 reinvent the > wheel? I hope not. I'm making some headway getting Lagrangespecific functions extended to handle more general elements; I'd rather not encourage the creation of more Lagrangeonly code!  Roy Stogner 
From: KIRK, BENJAMIN (JSCEG) (NASA) <benjamin.kirk1@na...>  20050321 15:20:01

It does now: MeshBase::element_iterator it = mesh.level_elements_begin(0); MeshBase::element_iterator end = mesh.level_elements_end(0); Ben Original Message From: libmeshusersadmin@... [mailto:libmeshusersadmin@...] On Behalf Of Michael Povolotskyi Sent: Monday, March 21, 2005 5:09 AM To: Roy Stogner Cc: libmeshusers@... Subject: Re: [Libmeshusers] Periodic boundary conditions > > 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.  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 _______________________________________________ Libmeshusers mailing list Libmeshusers@... https://lists.sourceforge.net/lists/listinfo/libmeshusers 
From: Michael Povolotskyi <povolotskyi@in...>  20050323 10:41:05

KIRK, BENJAMIN (JSCEG) (NASA) wrote: >It does now: > >MeshBase::element_iterator it = mesh.level_elements_begin(0); >MeshBase::element_iterator end = mesh.level_elements_end(0); > >Ben > > This code can not be compiled. I had to change to MeshBase::const_element_iterator el3 = mesh.level_elements_begin(0); MeshBase::const_element_iterator end_el3 = mesh.level_elements_end(0); Is this a correct solution? Michael. 
From: Benjamin S. Kirk <benkirk@cf...>  20050323 13:16:54

That is perfectly acceptable. I'm guessing in this case you have a const Mesh& and it won't let you access the nonconst build member. Ben Michael Povolotskyi wrote: > > > > KIRK, BENJAMIN (JSCEG) (NASA) wrote: > >> It does now: >> >> MeshBase::element_iterator it = mesh.level_elements_begin(0); >> MeshBase::element_iterator end = mesh.level_elements_end(0); >> >> Ben >> >> > This code can not be compiled. > I had to change to > > MeshBase::const_element_iterator el3 = mesh.level_elements_begin(0); > MeshBase::const_element_iterator end_el3 = mesh.level_elements_end(0); > > > Is this a correct solution? > Michael. > > > >  > This SF.net email is sponsored by: 2005 Windows Mobile Application Contest > Submit applications for Windows Mobile(tm)based Pocket PCs or Smartphones > for the chance to win $25,000 and application distribution. Enter today at > http://ads.osdn.com/?ad_id=6882&alloc_id=15148&op=click > _______________________________________________ > Libmeshusers mailing list > Libmeshusers@... > https://lists.sourceforge.net/lists/listinfo/libmeshusers 
From: John Peterson <peterson@cf...>  20050323 16:08:17

Michael Povolotskyi writes: > > > This code can not be compiled. > I had to change to > > MeshBase::const_element_iterator el3 = mesh.level_elements_begin(0); > MeshBase::const_element_iterator end_el3 = mesh.level_elements_end(0); > > > Is this a correct solution? If "mesh" is a constant reference in this case, then I believe that is correct. Consider #include <vector> void f(const std::vector<int>& v) { // std::vector<int>::iterator it = v.begin(); // error! std::vector<int>::const_iterator cit = v.begin(); // ok. } int main() { std::vector<int> v(10); f(v); return 0; } John 
From: Roy Stogner <roystgnr@ic...>  20100408 17:18:00

On Thu, 8 Apr 2010, David Knezevic wrote: >> Oh, and even if adaptivity is working, you'll definitely need the >> coarse mesh to match up. Otherwise there's no way to enforce >> periodic continuity without "locking" destroying your accuracy. > > Thanks for the info Roy. We're not doing adaptivity, we just generated a mesh > (using gmsh) for a periodic channel flow, but the mesh generator doesn't > necessarily match up inflow and outflow nodes so it sounds like that is a > problem? Definitely. Imagine a grid of piecewise linears or bilinears, intended to be periodic from left to right, and offset the nodes on the right side down by just a little bit. For strict C0 continuity, then, the slope on the left bottom element must match the slopes on the two right bottom elements. But the slope on the secondfrombottom element on the right must match the secondfrombottom on the left, etc., and by the time you're done your "piecewise bilinear" mesh has been constrained to a single linear along that entire edge. Disaster. Your only hope in that case is to enforce periodicity weakly: use a DG term or mortar method or some such on those boundaries, remembering that you'd have to add the appropriate terms to the send_list & sparsity pattern by hand. It would probably be easier to fix up the mesh generator output. Copying this to libmeshusers; I haven't documented the PeriodicBoundary usage well enough, and it would be nice if people Googling it found *something*...  Roy 
From: David Knezevic <dknez@MIT.EDU>  20100408 17:25:43

Roy Stogner wrote: > > > On Thu, 8 Apr 2010, David Knezevic wrote: > >>> Oh, and even if adaptivity is working, you'll definitely need the >>> coarse mesh to match up. Otherwise there's no way to enforce >>> periodic continuity without "locking" destroying your accuracy. >> >> Thanks for the info Roy. We're not doing adaptivity, we just generated >> a mesh (using gmsh) for a periodic channel flow, but the mesh >> generator doesn't necessarily match up inflow and outflow nodes so it >> sounds like that is a problem? > > Definitely. Imagine a grid of piecewise linears or bilinears, > intended to be periodic from left to right, and offset the nodes on > the right side down by just a little bit. For strict C0 continuity, > then, the slope on the left bottom element must match the slopes on > the two right bottom elements. But the slope on the > secondfrombottom element on the right must match the > secondfrombottom on the left, etc., and by the time you're done your > "piecewise bilinear" mesh has been constrained to a single linear > along that entire edge. Disaster. Very clear explanation, thanks. > Your only hope in that case is to enforce periodicity weakly: use a DG > term or mortar method or some such on those boundaries, remembering > that you'd have to add the appropriate terms to the send_list & > sparsity pattern by hand. It would probably be easier to fix up the > mesh generator output. Yes, the latter option is definitely easier. Dave 
From: Michael Povolotskyi <povolotskyi@in...>  20050315 14:34:27

Thank you, I will try this in a couple of days (several things has to be done before) and let you know how this works. In fact, I will need something more complicated because I want to use the adaptive mesh refinement, so the nodes at the opposite sides may not match. But I think that using interpolation similar to what is done for the hanging nodes I'll menage to impose periodic boundary conditions even in this case. Michael. KIRK, BENJAMIN (JSCEG) (NASA) wrote: >If you know the degree of freedom indices it should be easy. Essentially >what you want to do is replace one set of degrees of freedom with another. >Assuming u_m corresponds to DOFs n_u_m and similarly u_n is n_u_n you should >be able to do something like this: > >DofConstraintRow constraint; >constraint[n_u_n] = 1; >dof_map.add_constraint_row (n_u_m, constraint); > >Essentially this imposes u_m = u_n. Please let me know how this works out. > >Ben > > > > > > >Original Message >From: libmeshusersadmin@... >[mailto:libmeshusersadmin@...] On Behalf Of Michael >Povolotskyi >Sent: Tuesday, March 15, 2005 6:52 AM >To: libmesh  ML >Subject: [Libmeshusers] Periodic boundary conditions > > >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. > > > > > >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 >_______________________________________________ >Libmeshusers mailing list >Libmeshusers@... >https://lists.sourceforge.net/lists/listinfo/libmeshusers > > >   Michael Povolotskyi, Ph.D. University of Rome "Tor Vergata" Department of Electronic Engineering Viale Politecnico, 1  00133 Rome  Italy Phone + 39 06 72597367 Fax + 39 06 2020519 http://www.optolab.uniroma2.it/pages/moshe/moshe.html  
From: Roy Stogner <roystgnr@ic...>  20050315 20:08:49

On Tue, 15 Mar 2005, Michael Povolotskyi wrote: > In fact, I will need something more complicated because I want to > use the adaptive mesh refinement, so the nodes at the opposite sides > may not match. But I think that using interpolation similar to what > is done for the hanging nodes I'll menage to impose periodic > boundary conditions even in this case. Can I make a suggestion? Add the periodic boundary condition via penalty method: add the term "1/epsilon * integral over left boundary (left side value  right side value)" to your residual, and the corresponding matrix to your Jacobian. That'll be the easiest to code, should work for any problem with solutions in H1, and will work for any level of adaptivity with any finite element type.  Roy Stogner 
From: Michael Povolotskyi <povolotskyi@in...>  20050316 11:24:13

Thank you for the suggestion, I'll try the approach. However, I guess that in order to calculate "1/epsilon * integral over left boundary (left side value  right side value)" it will be necessary to do the same interpolation as has to be done to build a constraint. Michael. Roy Stogner wrote: > On Tue, 15 Mar 2005, Michael Povolotskyi wrote: > >> In fact, I will need something more complicated because I want to >> use the adaptive mesh refinement, so the nodes at the opposite sides >> may not match. But I think that using interpolation similar to what >> is done for the hanging nodes I'll menage to impose periodic >> boundary conditions even in this case. > > > Can I make a suggestion? Add the periodic boundary condition via > penalty method: add the term "1/epsilon * integral over left boundary > (left side value  right side value)" to your residual, and the > corresponding matrix to your Jacobian. That'll be the easiest to > code, should work for any problem with solutions in H1, and will work > for any level of adaptivity with any finite element type. >  > Roy Stogner > 
From: Roy Stogner <roystgnr@ic...>  20050316 13:33:10

On Wed, 16 Mar 2005, Michael Povolotskyi wrote: > Thank you for the suggestion, > I'll try the approach. However, I guess that in order to calculate > "1/epsilon * integral over left boundary > (left side value  right side value)" it will be necessary to do the same > interpolation as has to be done to build a constraint. Almost. To build a constraint, then for each DoF on one side of the boundary you'll have to: map its associated point onto the other side, find which cell contains that point, and find which DoFs in that cell constrain it. With the penalty method you'll just have to do the first two steps.  Roy Stogner 