From: Roy S. <roy...@ic...> - 2007-05-07 17:29:05
|
I hope you don't mind me copying these back to the list - if we've got a real bug here, it needs more attention than I can give it alone. In fact, I think I'll Cc: to Ben's email and make sure everyone's paying attention - I think Ben wrote build_constrained_matrix. ;-) --- Roy ---------- Forwarded message ---------- Date: Sun, 6 May 2007 00:27:54 -0500 From: Shun Wang <shu...@ui...> To: Roy Stogner <roy...@ic...> Subject: Re: [Libmesh-users] AMR in 3D I tested DofMap::build_constraint_matrix(). for both the smallest elements with node A. In one element, node A is interpolated as .75*C+.25*D. But in the other element, node A is interpolated as 0.5*C+0.5*B. The result from the second interpolation overwrites the first one and leaves the wrong solution on A. I don't quite understand the option called_recursively in DofMap::build_constrained_matrix(C,elem_dofs,called_recursively). It seems that DofMap::enforce_constraints_exactly() and DofMap::constrain_element_matrix_and_vector() use the default value "false" for called_recursively. I'm not sure if I should changed it to "true" in my case if I have this kind of mesh. My code is very complicated at this moment. I don't think it's easy to give you the test case. But if you could give me some suggestions or so, that would be very helpful. Thank you! On 5/6/07, Roy Stogner <roy...@ic...> wrote: > > On Sat, 5 May 2007, Shun Wang wrote: > >> If you checked the attached the picture. > > The picture didn't come through, but I take it you're talking about a > situation where the mesh is level-1 regular on faces but there is a > level-2 irregularity on edges? > >> This is a local mesh I get from libMesh AMR. There, both A and B >> are hanging nodes, while A depends on B. So when I do both > dofmap.enforce_constraints_exactly(), A is interpolated from B and >> C, while B as a hanging node has value 0. Therefore, A gets a >> wrong value. > > Hmm... the bug I found Friday should only have affected meshes which > have level-2 irregularities on sides, not just edges. > > If A is interpolated from B and C, while B is interpolated from > (suppose) D and E, then enforce_constraints_exactly() is supposed to > recursively build up a constraint matrix which directly interpolates A > from C, D, and E. It's of course possible that there are bugs in this > process. > >> I am guessing that libMesh level one condition is enforced only on >> surface, i.e. the level difference of two elements that share surfaces > can >> not exceed one. > >> Is my conjecture right? How can I solve this problem? > > There are now options for enforcing level one regularity between edge > and/or node neighbors in addition to face neighbors. If level-2 edge > irregularities are causing problems for your particular formulation, > using either of those options should fix it. > > However, even if your mesh has *no* regularity, libMesh is still > supposed to ensure that your constrained finite element spaces have > the appropriate level of continuity. If you're seeing an exception to > that rule, it's a bug and we'd appreciate if you can give us more > details or a reproducable test case. > --- > Roy > > ------------------------------------------------------------------------- > This SF.net email is sponsored by DB2 Express > Download DB2 Express C - the FREE version of DB2 express and take > control of your XML. No limits. Just data. Click to get it now. > http://sourceforge.net/powerbar/db2/ > _______________________________________________ > Libmesh-users mailing list > Lib...@li... > https://lists.sourceforge.net/lists/listinfo/libmesh-users > |
From: Roy S. <roy...@ic...> - 2007-05-07 17:29:12
|
---------- Forwarded message ---------- Date: Sun, 6 May 2007 00:31:12 -0500 From: Shun Wang <shu...@ui...> To: Roy Stogner <roy...@ic...> Subject: Re: [Libmesh-users] AMR in 3D And how do I check out this new piece of code for level-one mesh refinement from cvs? I don't want to check out all the code, because I added some code somewhere else in libMesh. Thanks! On 5/6/07, Shun Wang <shu...@ui...> wrote: > > I tested DofMap::build_constraint_matrix(). for both the smallest elements > with node A. > In one element, node A is interpolated as .75*C+.25*D. But in the other > element, node A is interpolated as 0.5*C+0.5*B. > The result from the second interpolation overwrites the first one and > leaves the wrong solution on A. > > I don't quite understand the option called_recursively in > DofMap::build_constrained_matrix(C,elem_dofs,called_recursively). It seems > that DofMap::enforce_constraints_exactly() and > DofMap::constrain_element_matrix_and_vector() use the default value "false" > for called_recursively. I'm not sure if I should changed it to "true" in my > case if I have this kind of mesh. > > My code is very complicated at this moment. I don't think it's easy to > give you the test case. But if you could give me some suggestions or so, > that would be very helpful. Thank you! > > On 5/6/07, Roy Stogner <roy...@ic...> wrote: >> >> On Sat, 5 May 2007, Shun Wang wrote: >> >>> If you checked the attached the picture. >> >> The picture didn't come through, but I take it you're talking about a >> situation where the mesh is level-1 regular on faces but there is a >> level-2 irregularity on edges? >> >>> This is a local mesh I get from libMesh AMR. There, both A and B >>> are hanging nodes, while A depends on B. So when I do both > > > > >> dofmap.enforce_constraints_exactly(), A is interpolated from B and >>> C, while B as a hanging node has value 0. Therefore, A gets a >>> wrong value. >> >> Hmm... the bug I found Friday should only have affected meshes which >> have level-2 irregularities on sides, not just edges. >> >> If A is interpolated from B and C, while B is interpolated from >> (suppose) D and E, then enforce_constraints_exactly() is supposed to >> recursively build up a constraint matrix which directly interpolates A >> from C, D, and E. It's of course possible that there are bugs in this >> process. >> >>> I am guessing that libMesh level one condition is enforced only on >>> surface, i.e. the level difference of two elements that share >> surfaces can >>> not exceed one. >> >>> Is my conjecture right? How can I solve this problem? >> >> There are now options for enforcing level one regularity between edge >> and/or node neighbors in addition to face neighbors. If level-2 edge >> irregularities are causing problems for your particular formulation, >> using either of those options should fix it. >> >> However, even if your mesh has *no* regularity, libMesh is still >> supposed to ensure that your constrained finite element spaces have >> the appropriate level of continuity. If you're seeing an exception to >> that rule, it's a bug and we'd appreciate if you can give us more >> details or a reproducable test case. >> --- >> Roy >> >> ------------------------------------------------------------------------- >> >> This SF.net email is sponsored by DB2 Express >> Download DB2 Express C - the FREE version of DB2 express and take >> control of your XML. No limits. Just data. Click to get it now. >> http://sourceforge.net/powerbar/db2/ >> _______________________________________________ >> Libmesh-users mailing list >> Lib...@li... >> https://lists.sourceforge.net/lists/listinfo/libmesh-users >> > > |
From: Roy S. <roy...@ic...> - 2007-05-07 17:39:01
|
> ---------- Forwarded message ---------- > Date: Sun, 6 May 2007 00:31:12 -0500 > From: Shun Wang <shu...@ui...> > To: Roy Stogner <roy...@ic...> > Subject: Re: [Libmesh-users] AMR in 3D > > And how do I check out this new piece of code for level-one mesh refinement > from cvs? I don't want to check out all the code, because I added some code > somewhere else in libMesh. Thanks! You can do a CVS update on individual files by just giving the filename at the end of the command; however I'd be very wary of trying to update libMesh piecemeal. You'd probably be better off doing a complete cvs update and then manually resolving any conflicts between the libMesh changes and your own changes. Also, if any of your code additions can be merged into the main tree without adverse effects, send us a patch! If you could also give us a small test case to use as an example program, there's even a chance that we won't break your patch with future libMesh changes. ;-) --- Roy |
From: John P. <pet...@cf...> - 2007-05-07 21:44:52
|
Hi, I don't have a huge amount of additional information to add to the thread, but I definitely think Shun Wang is on to something... I ran a transient, linear heat conduction problem on two different meshes www.cfdlab.ae.utexas.edu/~peterson/mismatched_hexes.png The mesh on the left has a maximum edge mismatch of 1 level while the one on the right (you might have to scroll over) has an edge mismatch level of 2. This is a linear problem, and so should converge in 1 'nonlinear' iteration. Here are the results for the first solve on the 1-level mismatched mesh: nonlinear |R| linear |R| 0.000624475 4.25804e-22 4.91167e-17 It's a symmetric problem and converges in 1 iteration, as expected. For the 2-level mismatch mesh, it's quite a different story: linear convergence of Newton's method (for a linear heat conduction problem!) nonlinear |R| linear |R| 0.000616188 3.60638e-22 6.75975e-06 1.27404e-24 1.88506e-07 3.13132e-25 5.23973e-09 1.29476e-26 1.44318e-10 3.53689e-28 3.93687e-12 1.37323e-30 I'm hopeful we can get this issue sorted out quickly, but I'm in the "don't know enough about build_constraint_matrix()" camp as well. -John Roy Stogner writes: > > > ---------- Forwarded message ---------- > > Date: Sun, 6 May 2007 00:31:12 -0500 > > From: Shun Wang <shu...@ui...> > > To: Roy Stogner <roy...@ic...> > > Subject: Re: [Libmesh-users] AMR in 3D > > > > And how do I check out this new piece of code for level-one mesh refinement > > from cvs? I don't want to check out all the code, because I added some code > > somewhere else in libMesh. Thanks! > > You can do a CVS update on individual files by just giving the > filename at the end of the command; however I'd be very wary of trying > to update libMesh piecemeal. You'd probably be better off doing a > complete cvs update and then manually resolving any conflicts between > the libMesh changes and your own changes. > > Also, if any of your code additions can be merged into the main tree > without adverse effects, send us a patch! If you could also give us a > small test case to use as an example program, there's even a chance > that we won't break your patch with future libMesh changes. ;-) > --- > Roy |
From: Lorenzo B. <bot...@gm...> - 2007-05-08 14:15:50
|
The fix for build_constraint_matrix() works... I need no more to enforce nodal mismatch and the solution doesn't blow up (I use max_h_level = 2). The function limit_level_mismatch_at_edge(...) included in the CVS doesn't work... I'm trying something else but I don't find a solution that works for all the type of elements (P1, P2, ecc...). I don't know if it is still useful, but probably it could be used with max_mismatch = 2. Thanks all for the help! Lorenzo |
From: John P. <pet...@cf...> - 2007-05-08 14:25:21
|
Lorenzo Botti writes: > The function limit_level_mismatch_at_edge(...) included in the CVS doesn't > work... Doesn't work as in 'crashes the program' or as in doesn't prevent the Newton iterations from diverging? -J |
From: Shun W. <shu...@ui...> - 2007-05-08 16:08:45
|
I'm glad that this quick fix works for you too. That means it is a real bug. However, I'm hoping for a more efficient solution. On 5/8/07, Lorenzo Botti <bot...@gm...> wrote: > > The fix for build_constraint_matrix() works... I need no more to enforce > nodal mismatch and the solution doesn't blow up (I use max_h_level = 2). > The function limit_level_mismatch_at_edge(...) included in the CVS > doesn't work... I'm trying something else but I don't find a solution that > works for all the type of elements (P1, P2, ecc...). I don't know if it is > still useful, but probably it could be used with max_mismatch = 2. > Thanks all for the help! > > Lorenzo > > ------------------------------------------------------------------------- > This SF.net email is sponsored by DB2 Express > Download DB2 Express C - the FREE version of DB2 express and take > control of your XML. No limits. Just data. Click to get it now. > http://sourceforge.net/powerbar/db2/ > _______________________________________________ > Libmesh-users mailing list > Lib...@li... > https://lists.sourceforge.net/lists/listinfo/libmesh-users > > |
From: Roy S. <roy...@ic...> - 2007-05-08 14:26:52
|
On Tue, 8 May 2007, Lorenzo Botti wrote: > The fix for build_constraint_matrix() works... I need no more to enforce > nodal mismatch and the solution doesn't blow up (I use max_h_level = 2). Good to hear. Now all we need to do is figure out a fix that won't break again at a level 3 edge mismatch... ;-) > The function limit_level_mismatch_at_edge(...) included in the CVS > doesn't work... You're right - the code for that new function doesn't work as is. I'll let you know when we've got a fixed version for you to try. --- Roy |
From: Lorenzo B. <bot...@gm...> - 2007-05-08 14:47:33
|
>You're right - the code for that new function doesn't work as is. >I'll let you know when we've got a fixed version for you to try. I wrote this code but I think it works only if the children (refined elements) "connected" to the vertices of the edge are both refining... bool MeshRefinement::limit_level_mismatch_at_edge (const unsigned int max_mismatch) { bool flags_changed = false; // Vector holding the maximum element level that touches a node. std::vector<unsigned char> max_level_at_node (_mesh.n_nodes(), 0); std::vector<unsigned char> max_p_level_at_node (_mesh.n_nodes(), 0); // Loop over all the active elements & fill the vector { MeshBase::element_iterator elem_it = _mesh.active_elements_begin(); const MeshBase::element_iterator elem_end = _mesh.active_elements_end(); for (; elem_it != elem_end; ++elem_it) { const Elem* elem = *elem_it; const unsigned char elem_level = elem->level() + ((elem->refinement_flag() == Elem::REFINE) ? 1 : 0); const unsigned char elem_p_level = elem->p_level() + ((elem->p_refinement_flag() == Elem::REFINE) ? 1 : 0); // Set the max_level at each node for (unsigned int n=0; n<elem->n_nodes(); n++) { const unsigned int node_number = elem->node(n); assert (node_number < max_level_at_node.size()); max_level_at_node[node_number] = std::max (max_level_at_node[node_number], elem_level); max_p_level_at_node[node_number] = std::max (max_p_level_at_node[node_number], elem_p_level); } } } // Now loop over the active elements and flag the elements // who violate the requested level mismatch { MeshBase::element_iterator elem_it = _mesh.active_elements_begin(); const MeshBase::element_iterator elem_end = _mesh.active_elements_end(); for (; elem_it != elem_end; ++elem_it) { Elem* elem = *elem_it; const unsigned int elem_level = elem->level(); const unsigned int elem_p_level = elem->p_level(); // Skip the element if it is already fully flagged if (elem->refinement_flag() == Elem::REFINE && elem->p_refinement_flag() == Elem::REFINE) continue; // Loop over the nodes, check for possible mismatch for (unsigned int n=0; n<elem->n_edges(); n++) { AutoPtr<Elem> edge = elem->build_edge(n); unsigned int node0 = edge->node(0); unsigned int node1 = edge->node(edge->n_nodes()-1); if ( (elem_level + max_mismatch) < max_level_at_node[node0] && (elem_level + max_mismatch) < max_level_at_node[node1] && elem->refinement_flag() != Elem::REFINE) { elem->set_refinement_flag (Elem::REFINE); flags_changed = true; } if ( (elem_p_level + max_mismatch) < max_p_level_at_node[node0] && (elem_p_level + max_mismatch) < max_p_level_at_node[node1] && elem->p_refinement_flag() != Elem::REFINE) { elem->set_p_refinement_flag (Elem::REFINE); flags_changed = true; } } } } return flags_changed; } 2007/5/8, Roy Stogner <roy...@ic...>: > > On Tue, 8 May 2007, Lorenzo Botti wrote: > > > The fix for build_constraint_matrix() works... I need no more to enforce > > nodal mismatch and the solution doesn't blow up (I use max_h_level = 2). > > Good to hear. Now all we need to do is figure out a fix that won't > break again at a level 3 edge mismatch... ;-) > > > The function limit_level_mismatch_at_edge(...) included in the CVS > > doesn't work... > > You're right - the code for that new function doesn't work as is. > I'll let you know when we've got a fixed version for you to try. > --- > Roy > |
From: John P. <pet...@cf...> - 2007-05-08 15:00:57
|
I don't think this is correct in general. It gives the mid-edge node for a quadratic edge... which might be why it started working for you? > unsigned int node1 = edge->node(edge->n_nodes()-1); Lorenzo Botti writes: > >You're right - the code for that new function doesn't work as is. > >I'll let you know when we've got a fixed version for you to try. > > I wrote this code but I think it works only if the children (refined > elements) "connected" to the vertices of the edge are both refining... > > bool MeshRefinement::limit_level_mismatch_at_edge (const unsigned int > max_mismatch) > { > bool flags_changed = false; > > > // Vector holding the maximum element level that touches a node. > std::vector<unsigned char> max_level_at_node (_mesh.n_nodes(), 0); > std::vector<unsigned char> max_p_level_at_node (_mesh.n_nodes(), 0); > > > // Loop over all the active elements & fill the vector > { > MeshBase::element_iterator elem_it = > _mesh.active_elements_begin(); > const MeshBase::element_iterator elem_end = _mesh.active_elements_end(); > > for (; elem_it != elem_end; ++elem_it) > { > const Elem* elem = *elem_it; > const unsigned char elem_level = > elem->level() + ((elem->refinement_flag() == Elem::REFINE) ? 1 : 0); > const unsigned char elem_p_level = > elem->p_level() + ((elem->p_refinement_flag() == Elem::REFINE) ? 1 : > 0); > > // Set the max_level at each node > for (unsigned int n=0; n<elem->n_nodes(); n++) > { > const unsigned int node_number = elem->node(n); > > assert (node_number < max_level_at_node.size()); > > max_level_at_node[node_number] = > std::max (max_level_at_node[node_number], elem_level); > max_p_level_at_node[node_number] = > std::max (max_p_level_at_node[node_number], elem_p_level); > } > } > } > > // Now loop over the active elements and flag the elements > // who violate the requested level mismatch > > { > MeshBase::element_iterator elem_it = > _mesh.active_elements_begin(); > const MeshBase::element_iterator elem_end = _mesh.active_elements_end(); > > for (; elem_it != elem_end; ++elem_it) > { > Elem* elem = *elem_it; > const unsigned int elem_level = elem->level(); > const unsigned int elem_p_level = elem->p_level(); > > // Skip the element if it is already fully flagged > if (elem->refinement_flag() == Elem::REFINE && > elem->p_refinement_flag() == Elem::REFINE) > continue; > > // Loop over the nodes, check for possible mismatch > for (unsigned int n=0; n<elem->n_edges(); n++) > { > AutoPtr<Elem> edge = elem->build_edge(n); > unsigned int node0 = edge->node(0); > unsigned int node1 = edge->node(edge->n_nodes()-1); > > if ( (elem_level + max_mismatch) < max_level_at_node[node0] > && (elem_level + max_mismatch) < max_level_at_node[node1] > && elem->refinement_flag() != Elem::REFINE) > { > elem->set_refinement_flag (Elem::REFINE); > flags_changed = true; > } > if ( (elem_p_level + max_mismatch) < max_p_level_at_node[node0] > && (elem_p_level + max_mismatch) < max_p_level_at_node[node1] > && elem->p_refinement_flag() != Elem::REFINE) > { > elem->set_p_refinement_flag (Elem::REFINE); > flags_changed = true; > } > } > } > > } > > return flags_changed; > } > > > > 2007/5/8, Roy Stogner <roy...@ic...>: > > > > On Tue, 8 May 2007, Lorenzo Botti wrote: > > > > > The fix for build_constraint_matrix() works... I need no more to enforce > > > nodal mismatch and the solution doesn't blow up (I use max_h_level = 2). > > > > Good to hear. Now all we need to do is figure out a fix that won't > > break again at a level 3 edge mismatch... ;-) > > > > > The function limit_level_mismatch_at_edge(...) included in the CVS > > > doesn't work... > > > > You're right - the code for that new function doesn't work as is. > > I'll let you know when we've got a fixed version for you to try. > > --- > > Roy > > > ------------------------------------------------------------------------- > This SF.net email is sponsored by DB2 Express > Download DB2 Express C - the FREE version of DB2 express and take > control of your XML. No limits. Just data. Click to get it now. > http://sourceforge.net/powerbar/db2/_______________________________________________ > Libmesh-users mailing list > Lib...@li... > https://lists.sourceforge.net/lists/listinfo/libmesh-users |
From: Roy S. <roy...@ic...> - 2007-05-07 17:29:21
|
---------- Forwarded message ---------- Date: Sun, 6 May 2007 03:15:02 -0500 From: Shun Wang <shu...@ui...> To: Roy Stogner <roy...@ic...> Subject: Re: [Libmesh-users] AMR in 3D I found the bug. In this special case, DofMap::build_constrain_matrix() would go wrong. a----b----c-------------d Let's look at the case in a simple line and consider the element with nodes b,c. When we call build_constrain_matrix() with elem_dofs={bc}, up to line 745 of dof_map_constraints.C (see doc page<http://libmesh.sourceforge.net/doxygen/classDofMap.php#c6f698a546466fbd8a8b5ca355d54013>), we have new_elem_dofs = {abcd}, because b depends on {ac}, and c depends on {ad}. And in the constrain matrix C, the coefficients on row b are 0.5 on column a and 0.5 on column c. Now build_constrain_matrix() is called recursively with elem_dofs={abcd}, then it won't expand any more and return an empty constrain matrix Cnew. Then, the constrain matrix C will have wrong row b. So, I changed the code as following to correct this bug: 1. remove the "if" condition on lines 699 and 700 2. put a "if (!called_recursively)" condition for the recursive call between lines 745 to 753. I think there is better ways to implement build_constrain_matrix() without recursive calls. One way is to do some pre-processing on DofConstraints to eliminate recursive dependency (hanging nodes depend on other hanging nodes). I'm not sure if this will make trouble for other things, but it will certainly make the frequently called function DofMap::build_constrain_matrix() more efficient. I hope this is helpful. On 5/6/07, Shun Wang <shu...@ui...> wrote: > > And how do I check out this new piece of code for level-one mesh > refinement from cvs? I don't want to check out all the code, because I added > some code somewhere else in libMesh. Thanks! > > On 5/6/07, Shun Wang <shu...@ui...> wrote: >> >> I tested DofMap::build_constraint_matrix(). for both the smallest >> elements with node A. >> In one element, node A is interpolated as .75*C+.25*D. But in the other >> element, node A is interpolated as 0.5*C+0.5*B. >> The result from the second interpolation overwrites the first one and >> leaves the wrong solution on A. >> >> I don't quite understand the option called_recursively in >> DofMap::build_constrained_matrix(C,elem_dofs,called_recursively). It seems >> that DofMap::enforce_constraints_exactly() and >> DofMap::constrain_element_matrix_and_vector() use the default value "false" >> for called_recursively. I'm not sure if I should changed it to "true" in my >> case if I have this kind of mesh. >> >> My code is very complicated at this moment. I don't think it's easy to >> give you the test case. But if you could give me some suggestions or so, >> that would be very helpful. Thank you! >> >> On 5/6/07, Roy Stogner <roy...@ic...> wrote: >>> >>> On Sat, 5 May 2007, Shun Wang wrote: >>> >>>> If you checked the attached the picture. >>> >>> The picture didn't come through, but I take it you're talking about a >>> situation where the mesh is level-1 regular on faces but there is a >>> level-2 irregularity on edges? >>> >>>> This is a local mesh I get from libMesh AMR. There, both A and B >>>> are hanging nodes, while A depends on B. So when I do both >> >> >> >> >>> dofmap.enforce_constraints_exactly(), A is interpolated from B and >>>> C, while B as a hanging node has value 0. Therefore, A gets a >>>> wrong value. >>> >>> Hmm... the bug I found Friday should only have affected meshes which >>> have level-2 irregularities on sides, not just edges. >>> >>> If A is interpolated from B and C, while B is interpolated from >>> (suppose) D and E, then enforce_constraints_exactly() is supposed to >>> recursively build up a constraint matrix which directly interpolates A >>> >>> from C, D, and E. It's of course possible that there are bugs in this >>> process. >>> >>>> I am guessing that libMesh level one condition is enforced only on >>>> surface, i.e. the level difference of two elements that share >>> surfaces can >>>> not exceed one. >>> >>>> Is my conjecture right? How can I solve this problem? >>> >>> There are now options for enforcing level one regularity between edge >>> and/or node neighbors in addition to face neighbors. If level-2 edge >>> irregularities are causing problems for your particular formulation, >>> using either of those options should fix it. >>> >>> However, even if your mesh has *no* regularity, libMesh is still >>> supposed to ensure that your constrained finite element spaces have >>> the appropriate level of continuity. If you're seeing an exception to >>> that rule, it's a bug and we'd appreciate if you can give us more >>> details or a reproducable test case. >>> --- >>> Roy >>> >>> ------------------------------------------------------------------------- >>> >>> This SF.net email is sponsored by DB2 Express >>> Download DB2 Express C - the FREE version of DB2 express and take >>> control of your XML. No limits. Just data. Click to get it now. >>> http://sourceforge.net/powerbar/db2/ >>> _______________________________________________ >>> Libmesh-users mailing list >>> Lib...@li... >>> https://lists.sourceforge.net/lists/listinfo/libmesh-users >>> >> >> > |
From: Roy S. <roy...@ic...> - 2007-05-07 18:14:18
|
> ---------- Forwarded message ---------- > Date: Sun, 6 May 2007 03:15:02 -0500 > From: Shun Wang <shu...@ui...> > To: Roy Stogner <roy...@ic...> > Subject: Re: [Libmesh-users] AMR in 3D > > I found the bug. In this special case, DofMap::build_constrain_matrix() > would go wrong. > > a----b----c-------------d > > Let's look at the case in a simple line and consider the element with nodes > b,c. > > When we call build_constrain_matrix() with elem_dofs={bc}, up to line 745 of > dof_map_constraints.C (see doc > page<http://libmesh.sourceforge.net/doxygen/classDofMap.php#c6f698a546466fbd8a8b5ca355d54013>), > we have > new_elem_dofs = {abcd}, because b depends on {ac}, and c depends on {ad}. > And in the constrain matrix C, the coefficients on row b are 0.5 on column a > and 0.5 on column c. > > Now build_constrain_matrix() is called recursively with elem_dofs={abcd}, > then it won't expand any more and return an empty constrain matrix Cnew. > Then, the constrain matrix C will have wrong row b. This analysis looks correct to me. Unless Ben has an explanation of something we're missing, I'd say you've found a real bug. > So, I changed the code as following to correct this bug: > 1. remove the "if" condition on lines 699 and 700 > 2. put a "if (!called_recursively)" condition for the recursive call between > lines 745 to 753. This fix doesn't look correct to me. With that if(!called_recursively) condition, you'd never get more than one level of recursion, and the code would work on level-2 mismatches but fail on level-3. > I think there is better ways to implement build_constrain_matrix() without > recursive calls. > One way is to do some pre-processing on DofConstraints to eliminate > recursive dependency (hanging nodes depend on other hanging nodes). I'm not > sure if this will make trouble for other things, but it will certainly make > the frequently called function DofMap::build_constrain_matrix() more > efficient. I like this idea, but it would be a little more tricky than a single pre-processing step. If you do the pr-eprocessing in DofMap::create_dof_constraints(), then whenever the user calls DofMap::add_constraint_row() there's a possibility that the constraint list will now have become recursive again. Fixing that by doing another set of pre-processing sweeps after every add_constraint_row() call would be too expensive, and I can't think of any other easy solutions. --- Roy |
From: Roy S. <roy...@ic...> - 2007-05-08 17:45:00
|
On Mon, 7 May 2007, Roy Stogner wrote: > I like this idea, but it would be a little more tricky than a single > pre-processing step. If you do the pr-eprocessing in > DofMap::create_dof_constraints(), then whenever the user calls > DofMap::add_constraint_row() there's a possibility that the constraint > list will now have become recursive again. Fixing that by doing > another set of pre-processing sweeps after every add_constraint_row() > call would be too expensive, and I can't think of any other easy > solutions. On closer examination, it looks like we're encouraging the user to add constraints only in a user_constrain() function, and in any case the constraints should all be added by the time System::init_data() gets called so that ImplicitSystem and EigenSystem classes can build sparsity patterns correctly. It seems like we could call a postprocessing step in System::init_data() and System::reinit() and be assured that no additional user constraints should be added afterward. --- Roy |
From: John P. <pet...@cf...> - 2007-05-08 19:21:24
|
Kirk, Benjamin (JSC-EG) writes: > I'm relying on y'all to help me make sense of these ideas and to post > them to the list. All of the sudden I cannot send anything to the > lists: > > < ndjsbar02.ndc.nasa.gov #5.0.0 X-Spam-Firewall; host > mail.sourceforge.net[66.35.250.206] said: 550-Postmaster verification > failed while checking <ben...@na...> 550-Called: > 198.120.25.85 550-Sent: RCPT TO:<pos...@na...> > 550-Response: 550 <pos...@na...>: Recipient address rejected: > local mail delivery is disabled 550-Several RFCs state that you are > required to have a postmaster 550-mailbox for each mail domain. This > host does not accept mail 550-from domains whose servers reject the > postmaster address. 550 Sender verify failed (in reply to RCPT TO > command)> > > > > In any case, after talking with John I think that the easiest fix for > all this nonsense it to remove the recursion out of the constraint > matrix construction. The idea I have is as follows: > > (1) build up the _dof_constraints in the standard way. As implemented, > this produces a whole series of level-1 constraints. > > (2) process the _dof_constraints such that there are no implicit > constraints. That is, if a=f(b,c) and b=g(c,d) > then perform the function composition such that a=f(g(c,d),c) direcely. > More clearly, modify all dof constraints as to be written in terms of > active DOFS only. > > (3) build the constraint matrix, there will be no need for recursion. > > > Now the complication gets shifted to step (2), but it should be much > more clear. Note that all the level-1 constraints must be constructed > first. Then I envision something like this > > for each row of _dof_constraints > for each entry in the row > if the constraining dof is itself constrained > replace the dof with its constraints scaled by the weight > end > until there are no constrained dofs in the row > end > > > > So return to the case that was mentioned yesterday: > a----b----c-------------d > > After step (1) > > _dof_constraints[b] = .5*a + .5*c > _dof_constraints[c] = .5*a + .5*d > > Now, in step 2, since c is constrained we need to modify the first row: > > _dof_constraints[b] = .5*a + .5*(.5*a + .5*d) > = .75*a + .25*d > > (and now we are done since neither a nor d are constrained.) > > > This should avoid the issue that was occuring before. The problem was > stopping the recursion when the size of an object did not change, when > in fact it should terminate when the constraint is formed only in terms > of active dofs. > > > Is this clear? I think implementing this as a step at the end of > DofMap::create_dof_constraints() should be fairly straightforward. > > Thoughts?? > > -Ben > > > > > > |
From: Michael S. <m-s...@us...> - 2007-05-08 22:04:20
|
Hello, it is fine that you are rewriting the constraint system. I remember to have done this quite some time ago because I happened to need user constraints and did not understand the recursive constraints. I propose to include inhomogeneous constraints. Up to now, only constraint equations of the type sum_i weight_i * dof_i = 0 are allowed. This will not even permit the implementation of a simple Dirichlet boundary condition using constraints. An inhomogeneous constraint would have a right-hand side here, sum_i weight_i * dof_i = f The way to implement this is, of course, to build the constraints for the matrix and the right-hand-side vector at the same time. Best Michael. On 08.05.07, John Peterson wrote: > Kirk, Benjamin (JSC-EG) writes: > > I'm relying on y'all to help me make sense of these ideas and to post > > them to the list. All of the sudden I cannot send anything to the > > lists: > > > > < ndjsbar02.ndc.nasa.gov #5.0.0 X-Spam-Firewall; host > > mail.sourceforge.net[66.35.250.206] said: 550-Postmaster verification > > failed while checking <ben...@na...> 550-Called: > > 198.120.25.85 550-Sent: RCPT TO:<pos...@na...> > > 550-Response: 550 <pos...@na...>: Recipient address rejected: > > local mail delivery is disabled 550-Several RFCs state that you are > > required to have a postmaster 550-mailbox for each mail domain. This > > host does not accept mail 550-from domains whose servers reject the > > postmaster address. 550 Sender verify failed (in reply to RCPT TO > > command)> > > > > > > > > In any case, after talking with John I think that the easiest fix for > > all this nonsense it to remove the recursion out of the constraint > > matrix construction. The idea I have is as follows: > > > > (1) build up the _dof_constraints in the standard way. As implemented, > > this produces a whole series of level-1 constraints. > > > > (2) process the _dof_constraints such that there are no implicit > > constraints. That is, if a=f(b,c) and b=g(c,d) > > then perform the function composition such that a=f(g(c,d),c) direcely. > > More clearly, modify all dof constraints as to be written in terms of > > active DOFS only. > > > > (3) build the constraint matrix, there will be no need for recursion. > > > > > > Now the complication gets shifted to step (2), but it should be much > > more clear. Note that all the level-1 constraints must be constructed > > first. Then I envision something like this > > > > for each row of _dof_constraints > > for each entry in the row > > if the constraining dof is itself constrained > > replace the dof with its constraints scaled by the weight > > end > > until there are no constrained dofs in the row > > end > > > > > > > > So return to the case that was mentioned yesterday: > > a----b----c-------------d > > > > After step (1) > > > > _dof_constraints[b] = .5*a + .5*c > > _dof_constraints[c] = .5*a + .5*d > > > > Now, in step 2, since c is constrained we need to modify the first row: > > > > _dof_constraints[b] = .5*a + .5*(.5*a + .5*d) > > = .75*a + .25*d > > > > (and now we are done since neither a nor d are constrained.) > > > > > > This should avoid the issue that was occuring before. The problem was > > stopping the recursion when the size of an object did not change, when > > in fact it should terminate when the constraint is formed only in terms > > of active dofs. > > > > > > Is this clear? I think implementing this as a step at the end of > > DofMap::create_dof_constraints() should be fairly straightforward. > > > > Thoughts?? > > > > -Ben > > > > > > > > > > > > > > ------------------------------------------------------------------------- > This SF.net email is sponsored by DB2 Express > Download DB2 Express C - the FREE version of DB2 express and take > control of your XML. No limits. Just data. Click to get it now. > http://sourceforge.net/powerbar/db2/ > _______________________________________________ > Libmesh-users mailing list > Lib...@li... > https://lists.sourceforge.net/lists/listinfo/libmesh-users -- Michael Schindler Laboratoire de Physico-Chimie Théorique. ESPCI. 10 rue Vauquelin, 75231 Paris cedex 05, France. Tel: +33 (0)1 40 79 58 41 Fax: +33 (0)1 40 79 47 31 http: www.pct.espci.fr/~michael |
From: Roy S. <roy...@ic...> - 2007-05-08 19:30:59
|
On Tue, 8 May 2007, John Peterson wrote: > Kirk, Benjamin (JSC-EG) writes: > > > This should avoid the issue that was occuring before. The problem was > > stopping the recursion when the size of an object did not change, when > > in fact it should terminate when the constraint is formed only in terms > > of active dofs. > > > > > > Is this clear? I think implementing this as a step at the end of > > DofMap::create_dof_constraints() should be fairly straightforward. We can't put it in create_dof_constraints(), or it'll break as soon as someone tries to add user constraints which either depend upon or are depended upon by hanging node constraints. We can call such a postprocessing function from init_data() and reinit(), though, and it should work. I'm working on the implementation now. --- Roy |
From: Roy S. <roy...@ic...> - 2007-05-08 22:24:49
|
On Tue, 8 May 2007, Roy Stogner wrote: > We can call such a postprocessing function from init_data() and > reinit(), though, and it should work. I'm working on the > implementation now. If anyone wants to try it out, the relevant files are attached. This gives the right debugging and solution output on the test cases I've tried, but I'd like to see some more positive feedback before committing it to CVS. --- Roy |
From: Roy S. <roy...@ic...> - 2007-05-08 22:33:53
|
On Tue, 8 May 2007, Roy Stogner wrote: > On Tue, 8 May 2007, Roy Stogner wrote: > >> We can call such a postprocessing function from init_data() and >> reinit(), though, and it should work. I'm working on the >> implementation now. > > If anyone wants to try it out, the relevant files are attached. This > gives the right debugging and solution output on the test cases I've > tried, but I'd like to see some more positive feedback before > committing it to CVS. The relevant files hit some kind of sourceforge filesize limit when attached; instead, you can find them here: http://cfdlab.ae.utexas.edu/~roystgnr/constraints/ --- Roy |
From: Shun W. <shu...@ui...> - 2007-05-09 04:45:57
|
Thanks. I'll try it tomorrow and get back to you. On 5/8/07, Roy Stogner <roy...@ic...> wrote: > > On Tue, 8 May 2007, Roy Stogner wrote: > > > On Tue, 8 May 2007, Roy Stogner wrote: > > > >> We can call such a postprocessing function from init_data() and > >> reinit(), though, and it should work. I'm working on the > >> implementation now. > > > > If anyone wants to try it out, the relevant files are attached. This > > gives the right debugging and solution output on the test cases I've > > tried, but I'd like to see some more positive feedback before > > committing it to CVS. > > The relevant files hit some kind of sourceforge filesize limit when > attached; instead, you can find them here: > > http://cfdlab.ae.utexas.edu/~roystgnr/constraints/ > --- > Roy > > ------------------------------------------------------------------------- > This SF.net email is sponsored by DB2 Express > Download DB2 Express C - the FREE version of DB2 express and take > control of your XML. No limits. Just data. Click to get it now. > http://sourceforge.net/powerbar/db2/ > _______________________________________________ > Libmesh-devel mailing list > Lib...@li... > https://lists.sourceforge.net/lists/listinfo/libmesh-devel > |
From: Lorenzo B. <bot...@gm...> - 2007-05-09 10:36:45
|
It works great... Lorenzo 2007/5/9, Roy Stogner <roy...@ic...>: > > On Tue, 8 May 2007, Roy Stogner wrote: > > > On Tue, 8 May 2007, Roy Stogner wrote: > > > >> We can call such a postprocessing function from init_data() and > >> reinit(), though, and it should work. I'm working on the > >> implementation now. > > > > If anyone wants to try it out, the relevant files are attached. This > > gives the right debugging and solution output on the test cases I've > > tried, but I'd like to see some more positive feedback before > > committing it to CVS. > > The relevant files hit some kind of sourceforge filesize limit when > attached; instead, you can find them here: > > http://cfdlab.ae.utexas.edu/~roystgnr/constraints/ > --- > Roy > > ------------------------------------------------------------------------- > This SF.net email is sponsored by DB2 Express > Download DB2 Express C - the FREE version of DB2 express and take > control of your XML. No limits. Just data. Click to get it now. > http://sourceforge.net/powerbar/db2/ > _______________________________________________ > Libmesh-users mailing list > Lib...@li... > https://lists.sourceforge.net/lists/listinfo/libmesh-users > |
From: Roy S. <roy...@ic...> - 2007-05-09 18:15:14
|
On Wed, 9 May 2007, Lorenzo Botti wrote: > It works great... No problems for John or I either, so I'm committing process_recursive_constraints() to CVS and turning it on. If anyone experiences problems let me know. We might also want to take the now-redundant recursion out of build_constraint_matrix() sometime, but it's shouldn't be hurting anything, so one code rewrite at a time. --- Roy |