From: Roy Stogner <roystgnr@ic...>  20070507 17:29:21

 Forwarded message  Date: Sun, 6 May 2007 03:15:02 0500 From: Shun Wang <shunwang@...> To: Roy Stogner <roystgnr@...> Subject: Re: [Libmeshusers] AMR in 3D I found the bug. In this special case, DofMap::build_constrain_matrix() would go wrong. abcd 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 preprocessing 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 <shunwang@...> wrote: > > And how do I check out this new piece of code for levelone 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 <shunwang@...> 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 <roystgnr@...> 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 level1 regular on faces but there is a >>> level2 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 level2 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 level2 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/ >>> _______________________________________________ >>> Libmeshusers mailing list >>> Libmeshusers@... >>> https://lists.sourceforge.net/lists/listinfo/libmeshusers >>> >> >> > 