You can subscribe to this list here.
2003 
_{Jan}

_{Feb}

_{Mar}

_{Apr}

_{May}

_{Jun}

_{Jul}

_{Aug}

_{Sep}
(2) 
_{Oct}
(2) 
_{Nov}
(27) 
_{Dec}
(31) 

2004 
_{Jan}
(6) 
_{Feb}
(15) 
_{Mar}
(33) 
_{Apr}
(10) 
_{May}
(46) 
_{Jun}
(11) 
_{Jul}
(21) 
_{Aug}
(15) 
_{Sep}
(13) 
_{Oct}
(23) 
_{Nov}
(1) 
_{Dec}
(8) 
2005 
_{Jan}
(27) 
_{Feb}
(57) 
_{Mar}
(86) 
_{Apr}
(23) 
_{May}
(37) 
_{Jun}
(34) 
_{Jul}
(24) 
_{Aug}
(17) 
_{Sep}
(50) 
_{Oct}
(24) 
_{Nov}
(10) 
_{Dec}
(60) 
2006 
_{Jan}
(47) 
_{Feb}
(46) 
_{Mar}
(127) 
_{Apr}
(19) 
_{May}
(26) 
_{Jun}
(62) 
_{Jul}
(47) 
_{Aug}
(51) 
_{Sep}
(61) 
_{Oct}
(42) 
_{Nov}
(50) 
_{Dec}
(33) 
2007 
_{Jan}
(60) 
_{Feb}
(55) 
_{Mar}
(77) 
_{Apr}
(102) 
_{May}
(82) 
_{Jun}
(102) 
_{Jul}
(169) 
_{Aug}
(117) 
_{Sep}
(80) 
_{Oct}
(37) 
_{Nov}
(51) 
_{Dec}
(43) 
2008 
_{Jan}
(71) 
_{Feb}
(94) 
_{Mar}
(98) 
_{Apr}
(125) 
_{May}
(54) 
_{Jun}
(119) 
_{Jul}
(60) 
_{Aug}
(111) 
_{Sep}
(118) 
_{Oct}
(125) 
_{Nov}
(119) 
_{Dec}
(94) 
2009 
_{Jan}
(109) 
_{Feb}
(38) 
_{Mar}
(93) 
_{Apr}
(88) 
_{May}
(29) 
_{Jun}
(57) 
_{Jul}
(53) 
_{Aug}
(48) 
_{Sep}
(68) 
_{Oct}
(151) 
_{Nov}
(23) 
_{Dec}
(35) 
2010 
_{Jan}
(84) 
_{Feb}
(60) 
_{Mar}
(184) 
_{Apr}
(112) 
_{May}
(60) 
_{Jun}
(90) 
_{Jul}
(23) 
_{Aug}
(70) 
_{Sep}
(119) 
_{Oct}
(27) 
_{Nov}
(47) 
_{Dec}
(54) 
2011 
_{Jan}
(22) 
_{Feb}
(19) 
_{Mar}
(92) 
_{Apr}
(93) 
_{May}
(35) 
_{Jun}
(91) 
_{Jul}
(32) 
_{Aug}
(61) 
_{Sep}
(7) 
_{Oct}
(69) 
_{Nov}
(81) 
_{Dec}
(23) 
2012 
_{Jan}
(64) 
_{Feb}
(95) 
_{Mar}
(35) 
_{Apr}
(36) 
_{May}
(63) 
_{Jun}
(98) 
_{Jul}
(70) 
_{Aug}
(171) 
_{Sep}
(149) 
_{Oct}
(64) 
_{Nov}
(67) 
_{Dec}
(126) 
2013 
_{Jan}
(108) 
_{Feb}
(104) 
_{Mar}
(171) 
_{Apr}
(133) 
_{May}
(108) 
_{Jun}
(100) 
_{Jul}
(93) 
_{Aug}
(126) 
_{Sep}
(74) 
_{Oct}
(59) 
_{Nov}
(145) 
_{Dec}
(93) 
2014 
_{Jan}
(38) 
_{Feb}
(45) 
_{Mar}
(26) 
_{Apr}
(41) 
_{May}
(125) 
_{Jun}
(70) 
_{Jul}
(61) 
_{Aug}
(64) 
_{Sep}

_{Oct}

_{Nov}

_{Dec}

S  M  T  W  T  F  S 



1
(6) 
2
(2) 
3
(1) 
4
(2) 
5

6

7
(3) 
8
(2) 
9
(2) 
10

11
(2) 
12

13

14
(7) 
15
(5) 
16
(1) 
17
(8) 
18
(2) 
19

20

21

22
(5) 
23
(5) 
24
(3) 
25
(1) 
26

27

28






From: Mark Blome <mblome@au...>  20050225 09:01:24

Hi everybody, I am currently working on a finite element application where I need to introduce additional degree of freedoms that are not associated to the mesh. I need them to approximate the region exterior to the finite element mesh. I already posted an email about this some weeks ago. Now I think I know (well more or less) how to create the matrix, in principle it looks something like:  Ke 0   x  =  Fe   0 K2   A   0  where Ke*x=Fe are the usual FEM equations and matrix K2 are additional entries and the vector A are the additional DOFs I need. Would it be a stupid idea to introduce something like a "large virtual finite element" that geometrically does not exist but has all the additional DOFs A associated to it? How could I do this? Is there any other more elegantly way of doing this in libmesh? I would be very happy about any comments on this, Mark 
From: KIRK, BENJAMIN (JSCEG) (NASA) <benjamin.kirk1@na...>  20050224 15:25:09

That is the correct behavior. As implemented, each system numbers its = DOFs independently into contiguous blocks. The reason for this is because, conceptually (at least thusfar) each system corresponds to a coupled = set of variables to be solved as a single linear system. Ben Original Message From: libmeshusersadmin@... [mailto:libmeshusersadmin@...] On Behalf Of Michael Schindler Sent: Thursday, February 24, 2005 4:22 AM To: libmeshusers@... Subject: [Libmeshusers] Problems with DofObject::dof_number Hello list, I have a couple of nodes and would like to know their dof numbers in = several systems: I have 3 systems, one with 3 variables, the others with one variable. Now, if I go through my list of (boundary) nodes and ask them = for their dofnumbers, I get the result, that they always report the = _first_ system, ignoring the other two. This is what I get by outputting the = system, variable, component and dof_number: system variable component > dof_number 0 0 0 > 1570 0 1 0 > 3143 0 2 0 > 3564 first node 1 0 0 > 1570 2 0 0 > 1570 0 0 0 > 1114 0 1 0 > 2687 0 2 0 > 3471 second node 1 0 0 > 1114 2 0 0 > 1114 0 0 0 > 1246 0 1 0 > 2819 0 2 0 > 3484 another node 1 0 0 > 1246 2 0 0 > 1246 This has been produced by the something like the following code: const Node * node =3D *nodes_itr; unsigned int n_sys =3D node>n_systems(); // all systems: loop begin for (sys_i=3D0; sys_i<n_sys; sys_i++) { n_vars =3D node>n_vars(sys_i); // all variables in this system: loop begin for (var_i=3D0; var_i<n_vars; var_i++) { n_comp =3D node>n_comp(sys_i, var_i); // components: loop begin for (comp_i=3D0; comp_i<n_comp; comp_i++) { dof_gi =3D node>dof_number(sys_i, var_i, comp_i); assert( dof_gi !=3D libMesh::invalid_uint ); std::cout << sys_i << " " << var_i << " " << comp_i << " > " = << cdof_gi << " " << dof_gi << std::endl; } // components: loop end } // variables: loop end } // systems: loop end I do not manipulate the dofs before this point (as far as I can see). Best greetings, Michael =20 "A mathematician is a device for turning coffee into theorems" Paul Erd=F6s.  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=3D6595&alloc_id=3D14396&op=3Dclick _______________________________________________ Libmeshusers mailing list Libmeshusers@... https://lists.sourceforge.net/lists/listinfo/libmeshusers 
From: Michael Schindler <mschindler@us...>  20050224 14:35:57

The problem has dissappeared  and showed to be no problem at all, but a stupid mistake. Of course, if I have several systems, the global dofs do not know of any other dof of a different system. Sorry for annoying you, Michael. On 24.02.05, Michael Schindler wrote: > Hello list, > > I have a couple of nodes and would like to know their dof numbers in > several systems: I have 3 systems, one with 3 variables, the others > with one variable. > Now, if I go through my list of (boundary) nodes and ask them for > their dofnumbers, I get the result, that they always report the > _first_ system, ignoring the other two. This is what I get by > outputting the system, variable, component and dof_number:  "A mathematician is a device for turning coffee into theorems" Paul Erdös. 
From: Michael Schindler <mschindler@us...>  20050224 10:22:36

Hello list, I have a couple of nodes and would like to know their dof numbers in several systems: I have 3 systems, one with 3 variables, the others with one variable. Now, if I go through my list of (boundary) nodes and ask them for their dofnumbers, I get the result, that they always report the _first_ system, ignoring the other two. This is what I get by outputting the system, variable, component and dof_number: system variable component > dof_number 0 0 0 > 1570 0 1 0 > 3143 0 2 0 > 3564 first node 1 0 0 > 1570 2 0 0 > 1570 0 0 0 > 1114 0 1 0 > 2687 0 2 0 > 3471 second node 1 0 0 > 1114 2 0 0 > 1114 0 0 0 > 1246 0 1 0 > 2819 0 2 0 > 3484 another node 1 0 0 > 1246 2 0 0 > 1246 This has been produced by the something like the following code: const Node * node = *nodes_itr; unsigned int n_sys = node>n_systems(); // all systems: loop begin for (sys_i=0; sys_i<n_sys; sys_i++) { n_vars = node>n_vars(sys_i); // all variables in this system: loop begin for (var_i=0; var_i<n_vars; var_i++) { n_comp = node>n_comp(sys_i, var_i); // components: loop begin for (comp_i=0; comp_i<n_comp; comp_i++) { dof_gi = node>dof_number(sys_i, var_i, comp_i); assert( dof_gi != libMesh::invalid_uint ); std::cout << sys_i << " " << var_i << " " << comp_i << " > " << cdof_gi << " " << dof_gi << std::endl; } // components: loop end } // variables: loop end } // systems: loop end I do not manipulate the dofs before this point (as far as I can see). Best greetings, Michael  "A mathematician is a device for turning coffee into theorems" Paul Erdös. 
From: Roy Stogner <roystgnr@ic...>  20050223 15:03:55

On Wed, 23 Feb 2005, Michael Schindler wrote: > I have also tried the MKL BLAS routines, together with a petsc 2.2.1 > linked to the MKL. This also works fine for the case that has failed > before. Ben's apparantly right: it's a limitation with the default Debian BLAS. http://bugs.debian.org/cgibin/bugreport.cgi?bug=283731 I hesitate to call it a BLAS bug; based on the results on that page it sounds like a problem in PETSc which other BLAS implementations are just less sensitive to.  Roy 
From: John Peterson <peterson@cf...>  20050223 14:47:55

Michael Schindler writes: > > The error does not occur when using a different solver with ksp_type > It only occured with (ksp_type gmres), not for solvers tfqmr, bcgs, > and bicg. > > I have also tried the MKL BLAS routines, together with a petsc 2.2.1 > linked to the MKL. This also works fine for the case that has failed > before. In that case, I vote for a bad version of GMRES in the PETSc/BLAS distribution for Debian systems. I would not venture to guess why it only shows up for some grids and not others. As an added bonus, you now have the most uptodate version of PETSc and some fancy BLAS installed :) John 
From: Michael Schindler <mschindler@us...>  20050223 14:38:39

Hello, many thanks for the quick answers. On 22.02.05, Roy Stogner wrote: > On Tue, 22 Feb 2005, John Peterson wrote: > > >> I get an error message: > >> > >> ** On entry to DGEMV parameter number 6 had an illegal value > > > >That's definitely an error from the BLAS, but one that I've never > >seen. > > It's definitely an error from the BLAS, but not necessarily an error > in the BLASusing code. I've seen Laspack segfaults that were merely > symptomatic of errors I had made editing DoFMap code; it's likely > something similar is occuring here. > > Michael, I assume you're running in debug (METHOD=dbg) mode? If not, > that should be the first thing you try  it's not guaranteed to catch > errors as soon as they occur, but it's usually better at it than > optimized mode. Indeed, I am using (METHOD=dbg) and the appropriate value for (BOPT=g_c++). The error does not occur when using a different solver with ksp_type It only occured with (ksp_type gmres), not for solvers tfqmr, bcgs, and bicg. I have also tried the MKL BLAS routines, together with a petsc 2.2.1 linked to the MKL. This also works fine for the case that has failed before. To summarize this, I am not really sure about where the error has come from. I did no changes in the DofMap. Maybe there is something fishy in the dofconstraints system I have written to replace the usual dof constraints. But I doubt this, because it uses only the standard public interface for the matrices, dofmap, ... Best greetings, Michael.  "A mathematician is a device for turning coffee into theorems" Paul Erdös. 
From: Michael Povolotskyi <povolotskyi@in...>  20050223 12:50:02

Hello everybody, we're going to solve nonlinear Poisson equation using the adaptive mesh refinement tools. We know that the classes for nonlinear solvers are not finished yet. However, we would like to try to solve something. For that we need to clarify some points related to the "hanging" node problem. In example #10 there is a method dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices). This method may reduce the size of the matrix Ke and Fe since only the independent degrees of freedom have to be left. After the linear system is solved one has to return back to DOF's with hanging nodes. Where and how is this done? Thank you, Michael. 
From: Michael Povolotskyi <povolotskyi@in...>  20050223 12:50:01

Yes, the fix works for a linear equation we want to solve. For the moment we have some problems with nonlinear implementation so it is difficult to say if the fix will be ok also for a nonlinear case. Thank you for the reference. Michael. KIRK, BENJAMIN (JSCEG) (NASA) wrote: > Imagine a refined edge in 2D normal to a Dirichletvalued surface. The > hanging node on this edge has its degrees of freedom constrained in terms of > the two adjacent test/trial functions, one of which lives on the Dirichlet > boundary. As such, the corresponding test function should be 0, but this is > not imposed in any way by the DofMap. > > I think I'll hold off on any implementation until after the FE Rodeo in 2 > weeks. It will be easier to hash through the details at a facetoface > meeting. I want to be sure that any changes will also work for the > nonLagrange case. Michael, can you use the penalty hack/fix for the time > being? > > As for a reference, see > "Computational Grids", G.F. Carey, Taylor & Francis, 1997. (specifically, > pg. 211) > > Perhaps I can write something up more formally and include it in the > documentation. > > Ben > > > > > > Original Message > From: Roy Stogner [mailto:roystgnr@...] > Sent: Thursday, February 17, 2005 7:32 PM > To: KIRK, BENJAMIN (JSCEG) (NASA) > Cc: libmeshusers@... > Subject: RE: [Libmeshusers] adaptive refinements with linear elements > > > On Thu, 17 Feb 2005, KIRK, BENJAMIN (JSCEG) (NASA) wrote: > > >>a nonrefined element the refined elements degreesoffreedom are >>constrained in terms of the nonrefined elements DOFs. What is >>happening in this case is the refined elements matrix contributions >>are being aliased back to the nonrefined element (as well they >>should). However, they never get constrained out as you would like >>them to because of the elementbyelement BC application process. > > > Should that be happening in a 2D code with Lagrange elements? It seems like > any constrained nodes in that case should be interior nodes. > > >>I have envisioned another set of DOF constraints that basically impose >>dirichlet values. I'll work on that over the weekend and hopefully >>get you something to try next week. > > > Would you describe the method you're implementing? I'd like something which > could be generalized to nonLagrange elements or fourth order problems if > possible... of course, I don't know if that is possible at all, but maybe > hearing your ideas would give me a few. >  > Roy >   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: KIRK, BENJAMIN (JSCEG) (NASA) <benjamin.kirk1@na...>  20050222 19:00:43

Imagine a refined edge in 2D normal to a Dirichletvalued surface. The hanging node on this edge has its degrees of freedom constrained in terms of the two adjacent test/trial functions, one of which lives on the Dirichlet boundary. As such, the corresponding test function should be 0, but this is not imposed in any way by the DofMap. I think I'll hold off on any implementation until after the FE Rodeo in 2 weeks. It will be easier to hash through the details at a facetoface meeting. I want to be sure that any changes will also work for the nonLagrange case. Michael, can you use the penalty hack/fix for the time being? As for a reference, see "Computational Grids", G.F. Carey, Taylor & Francis, 1997. (specifically, pg. 211) Perhaps I can write something up more formally and include it in the documentation. Ben Original Message From: Roy Stogner [mailto:roystgnr@...] Sent: Thursday, February 17, 2005 7:32 PM To: KIRK, BENJAMIN (JSCEG) (NASA) Cc: libmeshusers@... Subject: RE: [Libmeshusers] adaptive refinements with linear elements On Thu, 17 Feb 2005, KIRK, BENJAMIN (JSCEG) (NASA) wrote: > a nonrefined element the refined elements degreesoffreedom are > constrained in terms of the nonrefined elements DOFs. What is > happening in this case is the refined elements matrix contributions > are being aliased back to the nonrefined element (as well they > should). However, they never get constrained out as you would like > them to because of the elementbyelement BC application process. Should that be happening in a 2D code with Lagrange elements? It seems like any constrained nodes in that case should be interior nodes. > I have envisioned another set of DOF constraints that basically impose > dirichlet values. I'll work on that over the weekend and hopefully > get you something to try next week. Would you describe the method you're implementing? I'd like something which could be generalized to nonLagrange elements or fourth order problems if possible... of course, I don't know if that is possible at all, but maybe hearing your ideas would give me a few.  Roy 
From: Roy Stogner <roystgnr@ic...>  20050222 16:01:52

On Tue, 22 Feb 2005, John Peterson wrote: > Michael Schindler writes: > > > I am having a problem with large systems. > > When using a square mesh consisting of 10x10 bins (in 2D) > > my program works fine. When turning to 100x100 bins > > I get an error message: > > > > ** On entry to DGEMV parameter number 6 had an illegal value > > That's definitely an error from the BLAS, but one that I've never > seen. It's definitely an error from the BLAS, but not necessarily an error in the BLASusing code. I've seen Laspack segfaults that were merely symptomatic of errors I had made editing DoFMap code; it's likely something similar is occuring here. Michael, I assume you're running in debug (METHOD=dbg) mode? If not, that should be the first thing you try  it's not guaranteed to catch errors as soon as they occur, but it's usually better at it than optimized mode.  Roy Stogner 
From: KIRK, BENJAMIN (JSCEG) (NASA) <benjamin.kirk1@na...>  20050222 15:58:47

I've certainly used larger meshes in 2D without incident. I'm guessing that the Debianprovided BLAS are restricting the size of the system used in the matrixvector product routine. As John suggests, you can build a new PETSc from source and link with the Intel MKL BLAS (or the ATLAS BLAS, which you can compile from source). In fact, you can probably just install a new BLAS and hack your PETSc installation so that it uses yours instead. Try 'make echo_ldflags' to see the libraries that libMesh is linking with. This will help you find your system BLAS library to replace. In the mean time, you might run the problem serially using the laspack iterative solvers to make sure there isn't an issue there. Run './foo uselaspack' to force your application to use the laspack solvers. Ben Original Message From: libmeshusersadmin@... [mailto:libmeshusersadmin@...] On Behalf Of John Peterson Sent: Tuesday, February 22, 2005 8:55 AM To: Michael Schindler Cc: libmeshusers@... Subject: [Libmeshusers] Solver problems for big systems Michael Schindler writes: > > Hello list, > > I am having a problem with large systems. > When using a square mesh consisting of 10x10 bins (in 2D) > my program works fine. When turning to 100x100 bins > I get an error message: > > > ** On entry to DGEMV parameter number 6 had an illegal value That's definitely an error from the BLAS, but one that I've never seen. We use the Intel MKL BLAS mostly, and have not had many problems. It would probably be difficult for me to reproduce your error message. A 100x100 mesh isn't *that* large, so I don't think that should be the problem, although it is strange that you don't get the same message on the 10x10. If I recall correctly, some other people have had problems with Debian+PETSc before too. If you can, you might want to try building Petsc 2.2.1 from source and linking against the MKL BLAS from the intel web page. In the meantime, what happens if you try to solve the problem with Laspack solvers? What happens if you change the PETSc solver you are using with a ksp_type flag on the command line? John  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: John Peterson <peterson@cf...>  20050222 14:55:08

Michael Schindler writes: > > Hello list, > > I am having a problem with large systems. > When using a square mesh consisting of 10x10 bins (in 2D) > my program works fine. When turning to 100x100 bins > I get an error message: > > > ** On entry to DGEMV parameter number 6 had an illegal value That's definitely an error from the BLAS, but one that I've never seen. We use the Intel MKL BLAS mostly, and have not had many problems. It would probably be difficult for me to reproduce your error message. A 100x100 mesh isn't *that* large, so I don't think that should be the problem, although it is strange that you don't get the same message on the 10x10. If I recall correctly, some other people have had problems with Debian+PETSc before too. If you can, you might want to try building Petsc 2.2.1 from source and linking against the MKL BLAS from the intel web page. In the meantime, what happens if you try to solve the problem with Laspack solvers? What happens if you change the PETSc solver you are using with a ksp_type flag on the command line? John 
From: Michael Schindler <mschindler@us...>  20050222 12:59:38

Hello list, I am having a problem with large systems. When using a square mesh consisting of 10x10 bins (in 2D) my program works fine. When turning to 100x100 bins I get an error message: ** On entry to DGEMV parameter number 6 had an illegal value I am using libmesh 0.4.3 rc2 together with petsc 2.2.0 (the debian testing package) and blas 3 (which comes as a debian dependency for petsc). DGEMV seems to be a blas routine for matrixvector multiplication. I suspect that petsc is calling this. Can someone help me with this problem? Best greetings, Michael.  "A mathematician is a device for turning coffee into theorems" Paul Erdös. 
From: Michael Povolotskyi <povolotskyi@in...>  20050218 16:36:15

Hello Ben, thank you for the hint. In fact, our final goal is to solve a nonlinear Poisson equation. That's why we are interested in a "nonpenalty" methods and we are looking forward to the new DOF constraints. We'd greatly appreciate if you give us a reference that explains the DOF constraints that are already implemented in Libmesh. Michael. KIRK, BENJAMIN (JSCEG) (NASA) wrote: > I effectively reimplemented the penalty method in your code by changing > > Fe(i) += ... > Ke(i,j) += ... > > to > > Fe(i) += 1.e30*... > Ke(i,j) += 1.e30*... > > in the Dirichlet BC section, and everything worked out fine. > > The problem is that when you have an adaptively refined element neighboring > a nonrefined element the refined elements degreesoffreedom are > constrained in terms of the nonrefined elements DOFs. What is happening in > this case is the refined elements matrix contributions are being aliased > back to the nonrefined element (as well they should). However, they never > get constrained out as you would like them to because of the > elementbyelement BC application process. > > By using the penalty approach you circumvent this issue because of some > quirks in floatingpoint arithmetic (1.e30 + 1.e0 = 1.e30). However, there > should be a more convenient way to impose boundary conditions without > resorting to this approach. > > I have envisioned another set of DOF constraints that basically impose > dirichlet values. I'll work on that over the weekend and hopefully get you > something to try next week. In the meantime I advocate the above approach. > > Thanks for pointing this out! > > Ben > > > > > > Original Message > From: libmeshusersadmin@... > [mailto:libmeshusersadmin@...] On Behalf Of Michael > Povolotskyi > Sent: Thursday, February 17, 2005 9:31 AM > To: John Peterson > Cc: libmeshusers@... > Subject: Re: [Libmeshusers] adaptive refinements with linear elements > > > Hi, > this is the code (Poisson.C) and the input file (Poisson.in). > > As for the possibility to set BC in a wrong way ... of course, everything is > possible but it worked well with a single grid and with uniform grid > refinement. Michael. > > John Peterson wrote: > >>Hi, >> >>If possible, you should post your entire code to the list. I'm sure >>one of us could find any problems with boundary conditions, etc. >>fairly quickly that way. I think you mentioned you were not using the >>penalty method as in the examples. Is it possible you are not setting >>the BC's correctly on the triangular elements? >> >>John >> >> >>Michael Povolotskyi writes: >> > Dear Libmesh developers, >> > with your help we are getting some progress. >> > >> > As you have suggested we tried linear elements (Tri3 and Quad4). >> > In this case the countour lines look continuous everywhere. >> > >> > However, we still see some differences between Tri3 and Quad4 elements. >> > The solution with Quad4 looks very nicely. >> > >> > The solution with Tri3 has some problems. >> > Namely, there are several nodes that lie at the boundary with the > > boundary condition u = 0.0 whose > >> > value depends on the adaptive grid refinement. >> > >> > Without refinement the value is equal to 0.0, as it should be. >> > Then we did 5 adaptive refinements and the values were: >> > 0.004 1e15 1e15 1e27 0.001 . >> > >> > The futher refinements lead to a constant value at these nodes of about > > 0.002. > >> > >> > The value of about 0.001 seems to be too much for the numerical error. >> > >> > Thank you very much, >> > we hope that our feedback is usefull. >> > >> > Michael. >> > > >   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...>  20050218 01:32:27

On Thu, 17 Feb 2005, KIRK, BENJAMIN (JSCEG) (NASA) wrote: > a nonrefined element the refined elements degreesoffreedom are > constrained in terms of the nonrefined elements DOFs. What is happening in > this case is the refined elements matrix contributions are being aliased > back to the nonrefined element (as well they should). However, they never > get constrained out as you would like them to because of the > elementbyelement BC application process. Should that be happening in a 2D code with Lagrange elements? It seems like any constrained nodes in that case should be interior nodes. > I have envisioned another set of DOF constraints that basically impose > dirichlet values. I'll work on that over the weekend and hopefully get you > something to try next week. Would you describe the method you're implementing? I'd like something which could be generalized to nonLagrange elements or fourth order problems if possible... of course, I don't know if that is possible at all, but maybe hearing your ideas would give me a few.  Roy 
From: KIRK, BENJAMIN (JSCEG) (NASA) <benjamin.kirk1@na...>  20050217 23:09:00

I effectively reimplemented the penalty method in your code by changing Fe(i) += ... Ke(i,j) += ... to Fe(i) += 1.e30*... Ke(i,j) += 1.e30*... in the Dirichlet BC section, and everything worked out fine. The problem is that when you have an adaptively refined element neighboring a nonrefined element the refined elements degreesoffreedom are constrained in terms of the nonrefined elements DOFs. What is happening in this case is the refined elements matrix contributions are being aliased back to the nonrefined element (as well they should). However, they never get constrained out as you would like them to because of the elementbyelement BC application process. By using the penalty approach you circumvent this issue because of some quirks in floatingpoint arithmetic (1.e30 + 1.e0 = 1.e30). However, there should be a more convenient way to impose boundary conditions without resorting to this approach. I have envisioned another set of DOF constraints that basically impose dirichlet values. I'll work on that over the weekend and hopefully get you something to try next week. In the meantime I advocate the above approach. Thanks for pointing this out! Ben Original Message From: libmeshusersadmin@... [mailto:libmeshusersadmin@...] On Behalf Of Michael Povolotskyi Sent: Thursday, February 17, 2005 9:31 AM To: John Peterson Cc: libmeshusers@... Subject: Re: [Libmeshusers] adaptive refinements with linear elements Hi, this is the code (Poisson.C) and the input file (Poisson.in). As for the possibility to set BC in a wrong way ... of course, everything is possible but it worked well with a single grid and with uniform grid refinement. Michael. John Peterson wrote: > Hi, > > If possible, you should post your entire code to the list. I'm sure > one of us could find any problems with boundary conditions, etc. > fairly quickly that way. I think you mentioned you were not using the > penalty method as in the examples. Is it possible you are not setting > the BC's correctly on the triangular elements? > > John > > > Michael Povolotskyi writes: > > Dear Libmesh developers, > > with your help we are getting some progress. > > > > As you have suggested we tried linear elements (Tri3 and Quad4). > > In this case the countour lines look continuous everywhere. > > > > However, we still see some differences between Tri3 and Quad4 elements. > > The solution with Quad4 looks very nicely. > > > > The solution with Tri3 has some problems. > > Namely, there are several nodes that lie at the boundary with the boundary condition u = 0.0 whose > > value depends on the adaptive grid refinement. > > > > Without refinement the value is equal to 0.0, as it should be. > > Then we did 5 adaptive refinements and the values were: > > 0.004 1e15 1e15 1e27 0.001 . > > > > The futher refinements lead to a constant value at these nodes of about 0.002. > > > > The value of about 0.001 seems to be too much for the numerical error. > > > > Thank you very much, > > we hope that our feedback is usefull. > > > > Michael. >   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: Michael Povolotskyi <povolotskyi@in...>  20050217 18:00:40

Hello, the situation gets clearer. I have just checked the matrix. It turnes that: a) if I don't call dof_map.constrain_element_matrix_and_vector(), then the matrix comes out as I intend; b) if I call dof_map.constrain_element_matrix_and_vector(), then the corresponding raw has additional elements. In case a) the boundary nodes are okay, but the solution is discontinuous. In case b) the solution is continuous, but the boundary conditions are disturbed. The question is: is it possible somehow to have both continuous solution and apply boundary conditions in our way (not with a penalty method)? Thank you for your help, Michael. Roy Stogner wrote: > On Thu, 17 Feb 2005, Michael Povolotskyi wrote: > >> we set boundary condition in such a way that a row in the matrix that >> correspond to >> the boundary node i looks like this: >> >> 0 .... 0 1 0 .... 0, >> >> where 1 stands at the ith position. >> >> We apply the boundary condition while assembling the matrix. >> Not all the Dirichlet boundary points are wrong  just a few. >> And not always the same. > > > Have you checked the assembled system matrix, to make sure those rows > came out the way you intended? Have you checked the solution vector, > to make sure the Dirichlet node errors are there and not just another > GMVrelated glitch? Are you sure the nodes end up exactly right with > uniform refinement? >  > Roy Stogner >   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: Michael Povolotskyi <povolotskyi@in...>  20050217 15:41:18

Hi, we set boundary condition in such a way that a row in the matrix that correspond to the boundary node i looks like this: 0 .... 0 1 0 .... 0, where 1 stands at the ith position. We apply the boundary condition while assembling the matrix. Not all the Dirichlet boundary points are wrong  just a few. And not always the same. Michael. Roy Stogner wrote: > On Thu, 17 Feb 2005, Michael Povolotskyi wrote: > >> Namely, there are several nodes that lie at the boundary with the >> boundary condition u = 0.0 whose value depends on the adaptive grid >> refinement. >> >> Without refinement the value is equal to 0.0, as it should be. >> Then we did 5 adaptive refinements and the values were: >> 0.004 1e15 1e15 1e27 0.001 . >> >> The futher refinements lead to a constant value at these nodes of >> about 0.002. >> >> The value of about 0.001 seems to be too much for the numerical error. > > > How are you setting the boundary conditions? Are you reinforcing them > at every solve? Are all your Dirichlet boundary nodes wrong, and if > so are they wrong in the same direction? >  > Roy Stogner >   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: Michael Povolotskyi <povolotskyi@in...>  20050217 15:32:02

Hi, this is the code (Poisson.C) and the input file (Poisson.in). As for the possibility to set BC in a wrong way ... of course, everything is possible but it worked well with a single grid and with uniform grid refinement. Michael. John Peterson wrote: > Hi, > > If possible, you should post your entire code to the list. > I'm sure one of us could find any problems with boundary > conditions, etc. fairly quickly that way. I think you > mentioned you were not using the penalty method as in > the examples. Is it possible you are not setting the > BC's correctly on the triangular elements? > > John > > > Michael Povolotskyi writes: > > Dear Libmesh developers, > > with your help we are getting some progress. > > > > As you have suggested we tried linear elements (Tri3 and Quad4). > > In this case the countour lines look continuous everywhere. > > > > However, we still see some differences between Tri3 and Quad4 elements. > > The solution with Quad4 looks very nicely. > > > > The solution with Tri3 has some problems. > > Namely, there are several nodes that lie at the boundary with the boundary condition u = 0.0 whose > > value depends on the adaptive grid refinement. > > > > Without refinement the value is equal to 0.0, as it should be. > > Then we did 5 adaptive refinements and the values were: > > 0.004 1e15 1e15 1e27 0.001 . > > > > The futher refinements lead to a constant value at these nodes of about 0.002. > > > > The value of about 0.001 seems to be too much for the numerical error. > > > > Thank you very much, > > we hope that our feedback is usefull. > > > > Michael. >   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...>  20050217 15:21:46

On Thu, 17 Feb 2005, Michael Povolotskyi wrote: > Namely, there are several nodes that lie at the boundary with the boundary > condition u = 0.0 whose value depends on the adaptive grid refinement. > > Without refinement the value is equal to 0.0, as it should be. > Then we did 5 adaptive refinements and the values were: > 0.004 1e15 1e15 1e27 0.001 . > > The futher refinements lead to a constant value at these nodes of about > 0.002. > > The value of about 0.001 seems to be too much for the numerical error. How are you setting the boundary conditions? Are you reinforcing them at every solve? Are all your Dirichlet boundary nodes wrong, and if so are they wrong in the same direction?  Roy Stogner 
From: John Peterson <peterson@cf...>  20050217 14:47:02

Hi, If possible, you should post your entire code to the list. I'm sure one of us could find any problems with boundary conditions, etc. fairly quickly that way. I think you mentioned you were not using the penalty method as in the examples. Is it possible you are not setting the BC's correctly on the triangular elements? John Michael Povolotskyi writes: > Dear Libmesh developers, > with your help we are getting some progress. > > As you have suggested we tried linear elements (Tri3 and Quad4). > In this case the countour lines look continuous everywhere. > > However, we still see some differences between Tri3 and Quad4 elements. > The solution with Quad4 looks very nicely. > > The solution with Tri3 has some problems. > Namely, there are several nodes that lie at the boundary with the boundary condition u = 0.0 whose > value depends on the adaptive grid refinement. > > Without refinement the value is equal to 0.0, as it should be. > Then we did 5 adaptive refinements and the values were: > 0.004 1e15 1e15 1e27 0.001 . > > The futher refinements lead to a constant value at these nodes of about 0.002. > > The value of about 0.001 seems to be too much for the numerical error. > > Thank you very much, > we hope that our feedback is usefull. > > Michael. 
From: Michael Povolotskyi <povolotskyi@in...>  20050217 11:21:29

Dear Libmesh developers, with your help we are getting some progress. As you have suggested we tried linear elements (Tri3 and Quad4). In this case the countour lines look continuous everywhere. However, we still see some differences between Tri3 and Quad4 elements. The solution with Quad4 looks very nicely. The solution with Tri3 has some problems. Namely, there are several nodes that lie at the boundary with the boundary condition u = 0.0 whose value depends on the adaptive grid refinement. Without refinement the value is equal to 0.0, as it should be. Then we did 5 adaptive refinements and the values were: 0.004 1e15 1e15 1e27 0.001 . The futher refinements lead to a constant value at these nodes of about 0.002. The value of about 0.001 seems to be too much for the numerical error. Thank you very much, we hope that our feedback is usefull. Michael. 
From: KIRK, BENJAMIN (JSCEG) (NASA) <benjamin.kirk1@na...>  20050217 02:33:49

I will look through this in detail as soon as I can. We are in the middle of a simulated mission right now, so I can't even breathe until Friday. Sorry, but I'm sure we can straighten this out in relatively short order. Ben Original Message From: Michael Povolotskyi [mailto:povolotskyi@...] Sent: Wednesday, February 16, 2005 11:18 AM To: KIRK, BENJAMIN (JSCEG) (NASA) Cc: 'Roy Stogner'; libmeshusers@... Subject: Re: [Libmeshusers] hanging nodes and triangular elements Today we performed a set of tests in order to understand where the problem is (if any). We are solving the equation \nabla^2 u = exp((x^2 + y^2)) on a square 1.0<= x <= 1.0, 1.0<= y <= 1.0 The boundary conditions are as follows: u = 1.0 if 1.0<=x<=0.8, y = 1.0 u = 1.0 if 0.8<=x<=1.0, y = 1.0 u = 0.0 if 0.4<=x<=0.4, y = 1.0 For the rest of the boundary we assume natural boundary conditions (von Neumann). In order to implement the Dirichlet boundary conditions we DON'T use the penalty function method shown in the examples. Instead we assign boundary values to the nodes that lie on the boundary. We found that the code produces a reasonable solution if: a) there is no mesh refinement; b) if there is uniform mesh refinement; We have tested both triangular and rectangular finite elements, both first and second approximation order. The "problems" start to appear when we use the adaptive mesh refinement. Detailed study shows that the same problems occur BOTH with triangular and rectangular elements (in contrast to what I informed yesterday). We send a solution visualised by GMV program in files example_fullview.pdf and example_magnified.pdf. One can see in Fig. example_magnified.pdf that the solution at some points (indicated by red circles) is not continuous. We tried different solvers. The solution changes very slightly, but the problem persists. In order to be more specific we attach our code (Poisson.C) and the input file Poisson.in that the program reads. We don't understand if we are using libmesh in a wrong way or if there is a problem in libmesh itself. Thank you for your time, Michael. KIRK, BENJAMIN (JSCEG) (NASA) wrote: > Are you using CG as the iterative solver? I just realized that the > constraints are inserted into the system matrix in such a way that the > solution vector automatically satisfies the constraints, but this is > not symmetric. > > Try running your problem when CG and then with GMRES or BiCG. If the > latter two work then the aforementioned problem is probably the source > of your difficulty. Let me know & I'll get a patch together. > > Ben > > > > Original Message > From: libmeshusersadmin@... > [mailto:libmeshusersadmin@...] On Behalf Of Roy > Stogner > Sent: Tuesday, February 15, 2005 12:18 PM > To: Michael Povolotskyi > Cc: libmeshusers@... > Subject: Re: [Libmeshusers] hanging nodes and triangular elements > > > On Tue, 15 Feb 2005, Michael Povolotskyi wrote: > > >>Yes, the code works correctly with uniformly refined triangular >>elements. >> >>Could you, please, tell how can I check if the nodes are in a correct >>order? > > > If you're getting good results with uniform refinement, there's likely > nothing wrong with your coarse mesh  the distorted triangle problem > was just the only thing I'd ever seen that would give you problems > with a triangular mesh but not a quadrilateral mesh. > > Could you specify exactly what problems you're getting with adaptive > triangles that you aren't seeing with uniform triangles or adaptive > quads? Is it a bad convergence rate, a crash, or something else? >  > Roy Stogner > > >  > 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...>  20050216 18:24:43

On Wed, 16 Feb 2005, Michael Povolotskyi wrote: > We have tested both triangular and rectangular finite elements, both first > and second approximation order. > > The "problems" start to appear when we use the adaptive mesh refinement. > Detailed study shows that the same problems occur BOTH with triangular and > rectangular elements (in contrast to what I informed yesterday). > > We send a solution visualised by GMV program in files example_fullview.pdf > and example_magnified.pdf. One can see in Fig. example_magnified.pdf that > the solution at some points (indicated by red circles) is not continuous. > We tried different solvers. The solution changes very slightly, but the > problem persists. Are you sure you see this problem with linear and bilinear elements? With quadratic and biquadratic elements, you may be encountering a known problem with GMV output of solutions on higher order elements. libMesh outputs the same linear or bilinear elements, regardless of what your solution degree is. For uniform meshes this isn't as obvious since the result is still continuous; for adaptive meshes, however, a hanging node's value will be plotted on fine element neighbors but not coarse element neighbors, and unless the hanging node's value is precisely the (bi)linear interpolant between the other nodes on its edge/face, the resulting plot will be discontinuous.  Roy 