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}
(66) 
_{Sep}
(60) 
_{Oct}
(110) 
_{Nov}
(27) 
_{Dec}
(30) 
2015 
_{Jan}
(43) 
_{Feb}
(67) 
_{Mar}
(71) 
_{Apr}
(92) 
_{May}
(39) 
_{Jun}
(15) 
_{Jul}
(46) 
_{Aug}
(63) 
_{Sep}
(84) 
_{Oct}
(82) 
_{Nov}
(69) 
_{Dec}
(45) 
2016 
_{Jan}
(92) 
_{Feb}
(91) 
_{Mar}
(148) 
_{Apr}
(43) 
_{May}
(58) 
_{Jun}
(117) 
_{Jul}
(87) 
_{Aug}

_{Sep}

_{Oct}

_{Nov}

_{Dec}

S  M  T  W  T  F  S 


1

2
(3) 
3
(1) 
4
(1) 
5
(6) 
6
(2) 
7
(2) 
8
(5) 
9

10
(1) 
11
(2) 
12

13
(1) 
14

15
(7) 
16
(11) 
17
(6) 
18
(1) 
19
(1) 
20

21

22

23
(1) 
24

25
(1) 
26
(2) 
27
(1) 
28

29
(1) 
30
(4) 
31




From: Roy Stogner <roystgnr@ic...>  20070116 22:54:00

On Tue, 16 Jan 2007, Manav Bhatia wrote: > Is there a consistent way of handling the issue of applying dirichlet > boundary conditions on a mesh that is being refined/coarsened? > For instance, if a portion of the boundary has dirichlet bc applied to it, > and gets refined, then the new nodes created inside that portion also should > have the same bc applied to it. One way to do this is to check for nodes on > that boundary and apply bc on them. Do you have a recommendation on this? > > Similarly, with the case of Neumann bcs...? I usually just have some if statements or switch statements that get applied to each quadrature point on a boundary side. If your boundary conditions are applied on regions that don't precisely break down onto element boundaries, then that's what you have to do, but it's not a great idea for robust code, because it means your boundary conditions have to be hardcoded in C++. For more sophisticated software, what you want to do is set boundary_id values on boundary sides of your coarse mesh, then in your assembly routine call boundary_id() on each boundary side you iterate over. The library will make sure that the ids remain correct through refinement and coarsening. See example 13 or example 18.  Roy 
From: Manav Bhatia <manav@u.washington.edu>  20070116 22:40:34

Thanks Roy, for your reply. I have additional questions related to the aspect of boundary conditions with AMR. Is there a consistent way of handling the issue of applying dirichlet boundary conditions on a mesh that is being refined/coarsened? For instance, if a portion of the boundary has dirichlet bc applied to it, and gets refined, then the new nodes created inside that portion also should have the same bc applied to it. One way to do this is to check for nodes on that boundary and apply bc on them. Do you have a recommendation on this? Similarly, with the case of Neumann bcs...? Thanks, Manav On Jan 16, 2007, at 2:12 PM, Roy Stogner wrote: > On Tue, 16 Jan 2007, Manav Bhatia wrote: > >> Kindly you help me with a hint about the basic concept of the >> MeshSmoother class? >> In my understanding, it does not affect the boundary nodes, but >> changes some of the domain interior node locations after a >> refinement/ >> coarsening. What is the need for this? > > Mesh generators, as well as meshdistorting formulations for > highdeformation mechanics problems, often produce meshes with > elements that are more distorted than optimal for a finite element > solver, or even too distorted to solve on at all. Running such a mesh > through a good smoother algorithm is intended to fix that. > > A mesh smoother may affect boundary as well as interior nodes, just so > long as the domain of the mesh remains relatively unchanged. In > practice what I've seen is that smoothers will allow nodes to "slide" > along flat boundaries, but will keep nodes on corners, edges, or > curves pinned in place. >  > Roy 
From: Derek Gaston <friedmud@gm...>  20070116 22:25:57

Hmmm.... I finally got my office setup at home... I should really put some time into getting my mesh smoother stuff in a working state and checked in... The problem is that WoW: Burning Crusade came out today ;) We'll see which one wins in the fight for my time.... Derek On 1/16/07, Roy Stogner <roystgnr@...> wrote: > On Tue, 16 Jan 2007, Manav Bhatia wrote: > > > Kindly you help me with a hint about the basic concept of the > > MeshSmoother class? > > In my understanding, it does not affect the boundary nodes, but > > changes some of the domain interior node locations after a refinement/ > > coarsening. What is the need for this? > > Mesh generators, as well as meshdistorting formulations for > highdeformation mechanics problems, often produce meshes with > elements that are more distorted than optimal for a finite element > solver, or even too distorted to solve on at all. Running such a mesh > through a good smoother algorithm is intended to fix that. > > A mesh smoother may affect boundary as well as interior nodes, just so > long as the domain of the mesh remains relatively unchanged. In > practice what I've seen is that smoothers will allow nodes to "slide" > along flat boundaries, but will keep nodes on corners, edges, or > curves pinned in place. >  > Roy > >  > Take Surveys. Earn Cash. Influence the Future of IT > Join SourceForge.net's Techsay panel and you'll get the chance to share your > opinions on IT & business topics through brief surveys  and earn cash > http://www.techsay.com/default.php?page=join.php&p=sourceforge&CID=DEVDEV > _______________________________________________ > Libmeshusers mailing list > Libmeshusers@... > https://lists.sourceforge.net/lists/listinfo/libmeshusers > 
From: Roy Stogner <roystgnr@ic...>  20070116 22:11:38

On Tue, 16 Jan 2007, Manav Bhatia wrote: > Kindly you help me with a hint about the basic concept of the > MeshSmoother class? > In my understanding, it does not affect the boundary nodes, but > changes some of the domain interior node locations after a refinement/ > coarsening. What is the need for this? Mesh generators, as well as meshdistorting formulations for highdeformation mechanics problems, often produce meshes with elements that are more distorted than optimal for a finite element solver, or even too distorted to solve on at all. Running such a mesh through a good smoother algorithm is intended to fix that. A mesh smoother may affect boundary as well as interior nodes, just so long as the domain of the mesh remains relatively unchanged. In practice what I've seen is that smoothers will allow nodes to "slide" along flat boundaries, but will keep nodes on corners, edges, or curves pinned in place.  Roy 
From: Manav Bhatia <manav@u.washington.edu>  20070116 22:04:32

Hi, Kindly you help me with a hint about the basic concept of the MeshSmoother class? In my understanding, it does not affect the boundary nodes, but changes some of the domain interior node locations after a refinement/ coarsening. What is the need for this? Thanks, Manav 
From: Shane Blackett <s.blackett@au...>  20070116 19:51:56

Hi, That environment variable INCLUDE is being used by the scripts but it doesn't understand the spaces. You don't want the Visual Studio includes anyway so just 'unset INCLUDE'. We have cygwin working (although I haven't tried with mpi or petsc), just a serial, static build so far. Including Petsc would be useful but we just haven't done it so far. It is also on my todo list to build for mingw, which just means writing code for a bit of unix specific stuff for threads and timing I think. Shane. 
From: Fowler U. Frederik <ancibz@te...>  20070116 18:50:28

From: Roy Stogner <roystgnr@ic...>  20070116 17:48:32

On Tue, 16 Jan 2007, Tim Kr=F6ger wrote: > It would be nice if you at least tested whether this phenomenon also ha= ppens=20 > with your installation, in order to exclude version conflicts. It breaks for me unless I set pc_factor_zeropivot 0, but that happens with a lot of libMesh problems, where our penalty methods confuse PETSc. Try that option and see how it works for you. That may not be the only thing you need to tweak, though. For one thing, you may also want to try adding more linear solves. 250 isn't a general rule, it's just something that worked well for that one problem and default mesh size.  Roy 
From: <tim@ce...>  20070116 13:23:26

Dear John, On Tue, 16 Jan 2007, John Peterson wrote: > Upon closer look at your code I finally noticed this part: > What does it do exactly? Well, it's the obstacle. I.e., in this case, a (large) square area is considered to be solid and fixed, i.e. u=0 in this area. The idea is that this acts as a noslip boundary condition at the boundary of the obstacle. Within the obstacle, I'm not interested in the NavierStokes quantities anyway. However, removing these cells from the grid might be problematic in my real application for two reasons: First, I have a second equation that has to be solved in the whole domain, including the obstacles. Second, the size, number, and shape of obstacles changes in time. > It can be quite tricky trying to > work with two different "penalty" values, in case they accidentily > end up in the same row. In my opinion, they can't happen to be in the same row in this case (except for extremly rough grids) since the obstacle is away from the boundary of the domain. Also, the behaviour of the program does not change if 1e7 is used for *all* penalties. However, if 1e10 is used for all penalties, then the computation with obstacle fails already on the 20x20 grid. Best Regards, Tim 
From: John Peterson <peterson@cf...>  20070116 08:32:15

Hi Tim, Upon closer look at your code I finally noticed this part: What does it do exactly=3F It can be quite tricky trying to work with two different "penalty" values, in case they accidentily end up in the same row. /* qqq begin */ if(fabs(elem>centroid()(0)0.5)<=3D0.4 && fabs(elem>centroid()(= 1)0.5)<=3D0.4) =09{ =09 const Real penalty =3D 1.e7; =09 for (unsigned int i=3D0; i<n=5Fu=5Fdofs; i++) =09 { =09 Kuu(i,i) +=3D penalty; =09 Kvv(i,i) +=3D penalty; =09 } =09} /* qqq end */ Tim Kr=F6ger writes: > Dear Roy and John, >=20 > thank you for your nearly simultaneous answer. (Did you use a=20 > stopwatch to get that right=3F) (: >=20 > Actually, constraining away the onedimensional null space did not=20= > cure the problem in my case. Please find attached a modified versio= n=20 > of my test program that includes John's code snippet. Again, all th= e=20 > (very little) places where the program differs from ex13 are marked=20= > with "qqq". The behaviour is the same as before, i.e. the program=20= > works well without obstacle (rough and fine grid), and it works well= =20 > with obstacle on the rough grid (20x20), but it crashes with obstacl= e=20 > on the fine grid (100x100). I pinned the pressure to 100.0 at node = 0,=20 > but pinning it to 0.0 instead does not change the problem. >=20 > It would be nice if you at least tested whether this phenomenon also= =20 > happens with your installation, in order to exclude version conflict= s. >=20 > Thank you very much in advance, >=20 > Tim > /* $Id: ex13.C,v 1.16 2006/10/23 18:46:06 jwpeterson Exp $ */ >=20 > /* The Next Great Finite Element Library. */ > /* Copyright (C) 2003 Benjamin S. Kirk */ >=20 > /* This library is free software; you can redistribute it and/or */ > /* modify it under the terms of the GNU Lesser General Public */ > /* License as published by the Free Software Foundation; either */ > /* version 2.1 of the License, or (at your option) any later version= . */ >=20 > /* This library is distributed in the hope that it will be useful, *= / > /* but WITHOUT ANY WARRANTY; without even the implied warranty of */= > /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU= */ > /* Lesser General Public License for more details. */ >=20 > /* You should have received a copy of the GNU Lesser General Public = */ > /* License along with this library; if not, write to the Free Softwa= re */ > /* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 021111= 307 USA */ >=20 > // <h1>Example 13  Unsteady NavierStokes Equations  Unsteady Non= linear Systems of Equations</h1> > // > // This example shows how a simple, unsteady, nonlinear system of e= quations > // can be solved in parallel. The system of equations are the fami= liar > // NavierStokes equations for lowspeed incompressible fluid flow.= This > // example introduces the concept of the inner nonlinear loop for e= ach > // timestep, and requires a good deal of linear algebra numbercrun= ching > // at each step. If you have the General Mesh Viewer (GMV) install= ed, > // the script movie.sh in this directory will also take appropriate= screen > // shots of each of the solution files in the time sequence. These= rgb files > // can then be animated with the "animate" utility of ImageMagick i= f it is > // installed on your system. On a PIII 1GHz machine in debug mode,= this > // example takes a little over a minute to run. If you would like = to see > // a more detailed time history, or compute more timesteps, that is= certainly > // possible by changing the n=5Ftimesteps and dt variables below. >=20 > // C++ include files that we need > #include <iostream> > #include <algorithm> > #include <math.h> >=20 > // Basic include file needed for the mesh functionality. > #include "libmesh.h" > #include "mesh.h" > #include "mesh=5Fgeneration.h" > #include "gmv=5Fio.h" > #include "equation=5Fsystems.h" > #include "fe.h" > #include "quadrature=5Fgauss.h" > #include "dof=5Fmap.h" > #include "sparse=5Fmatrix.h" > #include "numeric=5Fvector.h" > #include "dense=5Fmatrix.h" > #include "dense=5Fvector.h" > #include "linear=5Fimplicit=5Fsystem.h" > #include "transient=5Fsystem.h" > #include "perf=5Flog.h" > #include "boundary=5Finfo.h" >=20 > // Some (older) compilers do not offer full stream=20 > // functionality, OStringStream works around this. > #include "o=5Fstring=5Fstream.h" >=20 > // For systems of equations the \p DenseSubMatrix > // and \p DenseSubVector provide convenient ways for > // assembling the element matrix and vector on a > // componentbycomponent basis. > #include "dense=5Fsubmatrix.h" > #include "dense=5Fsubvector.h" >=20 > // The definition of a geometric element > #include "elem.h" >=20 > // Function prototype. This function will assemble the system > // matrix and righthandside. > void assemble=5Fstokes (EquationSystems& es, > =09=09 const std::string& system=5Fname); >=20 > // The main program. > int main (int argc, char** argv) > { > // Initialize libMesh. > libMesh::init (argc, argv); > { =20 > // Set the dimensionality of the mesh =3D 2 > const unsigned int dim =3D 2; =20 > =20 > // Create a twodimensional mesh. > Mesh mesh (dim); > =20 > // Use the MeshTools::Generation mesh generator to create a unif= orm > // grid on the square [1,1]^D. We instruct the mesh generator > // to build a mesh of 8x8 \p Quad9 elements in 2D, or \p Hex27 > // elements in 3D. Building these higherorder elements allows > // us to use higherorder approximation, as in example 3. > MeshTools::Generation::build=5Fsquare (mesh, > =09=09=09=09=09 //qqq=09=09=09=09=09 20, 20, > =09=09=09=09=09 100, 100,//qqq > =09=09=09=09=09 0., 1., > =09=09=09=09=09 0., 1., > =09=09=09=09=09 QUAD9); > =20 > // Print information about the mesh to the screen. > mesh.print=5Finfo(); > =20 > // Create an equation systems object. > EquationSystems equation=5Fsystems (mesh); > =20 > // Declare the system and its variables. > { > // Creates a transient system named "NavierStokes" > TransientLinearImplicitSystem & system =3D=20 > =09equation=5Fsystems.add=5Fsystem<TransientLinearImplicitSystem> ("= NavierStokes"); > =20 > // Add the variables "u" & "v" to "NavierStokes". They > // will be approximated using secondorder approximation. > system.add=5Fvariable ("u", SECOND); > system.add=5Fvariable ("v", SECOND); >=20 > // Add the variable "p" to "NavierStokes". This will > // be approximated with a firstorder basis, > // providing an LBBstable pressurevelocity pair. > system.add=5Fvariable ("p", FIRST); >=20 > // Give the system a pointer to the matrix assembly > // function. > system.attach=5Fassemble=5Ffunction (assemble=5Fstokes); > =20 > // Initialize the data structures for the equation system. > equation=5Fsystems.init (); >=20 > equation=5Fsystems.parameters.set<unsigned int>("linear solver= maximum iterations") =3D 250; > equation=5Fsystems.parameters.set<Real> ("linear solver= tolerance") =3D std::sqrt(TOLERANCE); > =20 > // Prints information about the system to the screen. > equation=5Fsystems.print=5Finfo(); > } >=20 > // Create a performancelogging object for this example > PerfLog perf=5Flog("Example 13"); > =20 > // Now we begin the timestep loop to compute the timeaccurate > // solution of the equations. > const Real dt =3D 0.005; > Real time =3D 0.0; > const unsigned int n=5Ftimesteps =3D 15; >=20 > // The number of steps and the stopping criterion are also requi= red > // for the nonlinear iterations. > const unsigned int n=5Fnonlinear=5Fsteps =3D 15; > const Real nonlinear=5Ftolerance =3D 1.e2; > =20 > // Tell the system of equations what the timestep is by using > // the set=5Fparameter function. The matrix assembly routine ca= n > // then reference this parameter. > equation=5Fsystems.parameters.set<Real> ("dt") =3D dt; >=20 > // Get a reference to the Stokes system to use later. > TransientLinearImplicitSystem& navier=5Fstokes=5Fsystem =3D > =09 equation=5Fsystems.get=5Fsystem<TransientLinearImplicitSystem>(= "NavierStokes"); >=20 > // The first thing to do is to get a copy of the solution at > // the current nonlinear iteration. This value will be used to > // determine if we can exit the nonlinear loop. > AutoPtr<NumericVector<Number> > > last=5Fnonlinear=5Fsoln (navier=5Fstokes=5Fsystem.solution>cl= one()); >=20 > for (unsigned int t=5Fstep=3D0; t=5Fstep<n=5Ftimesteps; ++t=5Fst= ep) > { > =09// Incremenet the time counter, set the time and the > =09// time step size as parameters in the EquationSystem. > =09time +=3D dt; >=20 > =09// Let the system of equations know the current time. > =09// This might be necessary for a timedependent forcing > =09// function for example. > =09equation=5Fsystems.parameters.set<Real> ("time") =3D time; >=20 > =09// A pretty update message > =09std::cout << " Solving time step " << t=5Fstep << ", time =3D " <= < time << std::endl; >=20 > =09// Now we need to update the solution vector from the > =09// previous time step. This is done directly through > =09// the reference to the Stokes system. > =09*navier=5Fstokes=5Fsystem.old=5Flocal=5Fsolution =3D *navier=5Fst= okes=5Fsystem.current=5Flocal=5Fsolution; >=20 > =09// Now we begin the nonlinear loop > =09for (unsigned int l=3D0; l<n=5Fnonlinear=5Fsteps; ++l) > =09 { > =09 // Update the nonlinear solution. > =09 last=5Fnonlinear=5Fsoln>zero(); > =09 last=5Fnonlinear=5Fsoln>add(*navier=5Fstokes=5Fsystem.soluti= on); > =09 =20 > =09 // Assemble & solve the linear system. > =09 perf=5Flog.start=5Fevent("linear solve"); > =09 equation=5Fsystems.get=5Fsystem("NavierStokes").solve(); > =09 perf=5Flog.stop=5Fevent("linear solve"); >=20 > =09 // Compute the difference between this solution and the last > =09 // nonlinear iterate. > =09 last=5Fnonlinear=5Fsoln>add (1., *navier=5Fstokes=5Fsystem.= solution); >=20 > =09 // Close the vector before computing its norm > =09 last=5Fnonlinear=5Fsoln>close(); >=20 > =09 // Compute the l2 norm of the difference > =09 const Real norm=5Fdelta =3D last=5Fnonlinear=5Fsoln>l2=5Fnor= m(); >=20 > =09 // Print out convergence information for the linear and > =09 // nonlinear iterations. > =09 std::cout << "Linear solver converged at step: " > =09=09 << navier=5Fstokes=5Fsystem.n=5Flinear=5Fiterations() > =09=09 << ", final residual: " > =09=09 << navier=5Fstokes=5Fsystem.final=5Flinear=5Fresidual() > =09=09 << " Nonlinear convergence: u  u=5Fold =3D " > =09=09 << norm=5Fdelta > =09=09 << std::endl; >=20 > =09 // Terminate the solution iteration if the difference between= > =09 // this iteration and the last is sufficiently small. > =09 if (norm=5Fdelta < nonlinear=5Ftolerance) > =09 { > =09=09std::cout << " Nonlinear solver converged at step " > =09=09=09 << l > =09=09=09 << std::endl; > =09=09break; > =09 } > =09 =20 > =09 } // end nonlinear loop > =09 > =09// Write out every nth timestep to file. > =09const unsigned int write=5Finterval =3D 1; > =09 > =09if ((t=5Fstep+1)%write=5Finterval =3D=3D 0) > =09 { > =09 OStringStream file=5Fname; >=20 > =09 // We write the file name in the gmv autoread format. > =09 file=5Fname << "out.gmv."; > =09 OSSRealzeroright(file=5Fname,3,0, t=5Fstep + 1); > =09 =20 > =09 GMVIO(mesh).write=5Fequation=5Fsystems (file=5Fname.str(), > =09=09=09=09=09=09equation=5Fsystems); > =09 } > } // end timestep loop. > } > =20 > // All done. =20 > return libMesh::close (); > } >=20 >=20 >=20 >=20 >=20 >=20 > // The matrix assembly function to be called at each time step to > // prepare for the linear solve. > void assemble=5Fstokes (EquationSystems& es, > =09=09 const std::string& system=5Fname) > { > // It is a good idea to make sure we are assembling > // the proper system. > assert (system=5Fname =3D=3D "NavierStokes"); > =20 > // Get a constant reference to the mesh object. > const Mesh& mesh =3D es.get=5Fmesh(); > =20 > // The dimension that we are running > const unsigned int dim =3D mesh.mesh=5Fdimension(); > =20 > // Get a reference to the Stokes system object. > TransientLinearImplicitSystem & navier=5Fstokes=5Fsystem =3D > es.get=5Fsystem<TransientLinearImplicitSystem> ("NavierStokes")= ; >=20 > // Numeric ids corresponding to each variable in the system > const unsigned int u=5Fvar =3D navier=5Fstokes=5Fsystem.variable=5F= number ("u"); > const unsigned int v=5Fvar =3D navier=5Fstokes=5Fsystem.variable=5F= number ("v"); > const unsigned int p=5Fvar =3D navier=5Fstokes=5Fsystem.variable=5F= number ("p"); > =20 > // Get the Finite Element type for "u". Note this will be > // the same as the type for "v". > FEType fe=5Fvel=5Ftype =3D navier=5Fstokes=5Fsystem.variable=5Ftyp= e(u=5Fvar); > =20 > // Get the Finite Element type for "p". > FEType fe=5Fpres=5Ftype =3D navier=5Fstokes=5Fsystem.variable=5Fty= pe(p=5Fvar); >=20 > // Build a Finite Element object of the specified type for > // the velocity variables. > AutoPtr<FEBase> fe=5Fvel (FEBase::build(dim, fe=5Fvel=5Ftype)); > =20 > // Build a Finite Element object of the specified type for > // the pressure variables. > AutoPtr<FEBase> fe=5Fpres (FEBase::build(dim, fe=5Fpres=5Ftype)); > =20 > // A Gauss quadrature rule for numerical integration. > // Let the \p FEType object decide what order rule is appropriate.= > QGauss qrule (dim, fe=5Fvel=5Ftype.default=5Fquadrature=5Forder())= ; >=20 > // Tell the finite element objects to use our quadrature rule. > fe=5Fvel>attach=5Fquadrature=5Frule (&qrule); > fe=5Fpres>attach=5Fquadrature=5Frule (&qrule); > =20 > // Here we define some references to cellspecific data that > // will be used to assemble the linear system. > // > // The element Jacobian * quadrature weight at each integration po= int. =20 > const std::vector<Real>& JxW =3D fe=5Fvel>get=5FJxW(); >=20 > // The element shape functions evaluated at the quadrature points.= > const std::vector<std::vector<Real> >& phi =3D fe=5Fvel>get=5Fphi= (); >=20 > // The element shape function gradients for the velocity > // variables evaluated at the quadrature points. > const std::vector<std::vector<RealGradient> >& dphi =3D fe=5Fvel>= get=5Fdphi(); >=20 > // The element shape functions for the pressure variable > // evaluated at the quadrature points. > const std::vector<std::vector<Real> >& psi =3D fe=5Fpres>get=5Fph= i(); >=20 > // The value of the linear shape function gradients at the quadrat= ure points > // const std::vector<std::vector<RealGradient> >& dpsi =3D fe=5Fpr= es>get=5Fdphi(); > =20 > // A reference to the \p DofMap object for this system. The \p Do= fMap > // object handles the index translation from node and element numb= ers > // to degree of freedom numbers. We will talk more about the \p D= ofMap > // in future examples. > const DofMap & dof=5Fmap =3D navier=5Fstokes=5Fsystem.get=5Fdof=5F= map(); >=20 > // Define data structures to contain the element matrix > // and righthandside vector contribution. Following > // basic finite element terminology we will denote these > // "Ke" and "Fe". > DenseMatrix<Number> Ke; > DenseVector<Number> Fe; >=20 > DenseSubMatrix<Number> > Kuu(Ke), Kuv(Ke), Kup(Ke), > Kvu(Ke), Kvv(Ke), Kvp(Ke), > Kpu(Ke), Kpv(Ke), Kpp(Ke); >=20 > DenseSubVector<Number> > Fu(Fe), > Fv(Fe), > Fp(Fe); >=20 > // 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=5Findices; > std::vector<unsigned int> dof=5Findices=5Fu; > std::vector<unsigned int> dof=5Findices=5Fv; > std::vector<unsigned int> dof=5Findices=5Fp; >=20 > // Find out what the timestep size parameter is from the system, a= nd > // the value of theta for the theta method. We use implicit Euler= (theta=3D1) > // for this simulation even though it is only firstorder accurate= in time. > // The reason for this decision is that the secondorder CrankNic= olson > // method is notoriously oscillatory for problems with discontinuo= us > // initial data such as the liddriven cavity. Therefore, > // we sacrifice accuracy in time for stability, but since the solu= tion > // reaches steady state relatively quickly we can afford to take s= mall > // timesteps. If you monitor the initial nonlinear residual for t= his > // simulation, you should see that it is monotonically decreasing = in time. > const Real dt =3D es.parameters.get<Real>("dt"); > // const Real time =3D es.parameters.get<Real>("time"); > const Real theta =3D 1.; > =20 > // Now we will loop over all the elements in the mesh that > // live on the local processor. We will compute the element > // matrix and righthandside contribution. Since the mesh > // will be refined we want to only consider the ACTIVE elements, > // hence we use a variant of the \p active=5Felem=5Fiterator. > // const=5Factive=5Flocal=5Felem=5Fiterator el (mesh.ele= ments=5Fbegin()); > // const const=5Factive=5Flocal=5Felem=5Fiterator end=5Fel (mesh.e= lements=5Fend()); >=20 > MeshBase::const=5Felement=5Fiterator el =3D mesh.active=5F= local=5Felements=5Fbegin(); > const MeshBase::const=5Felement=5Fiterator end=5Fel =3D mesh.activ= e=5Flocal=5Felements=5Fend();=20 > =20 > for ( ; el !=3D end=5Fel; ++el) > { =20 > // Store a pointer to the element we are currently > // working on. This allows for nicer syntax later. > const Elem* elem =3D *el; > =20 > // Get the degree of freedom indices for the > // current element. These define where in the global > // matrix and righthandside this element will > // contribute to. > dof=5Fmap.dof=5Findices (elem, dof=5Findices); > dof=5Fmap.dof=5Findices (elem, dof=5Findices=5Fu, u=5Fvar); > dof=5Fmap.dof=5Findices (elem, dof=5Findices=5Fv, v=5Fvar); > dof=5Fmap.dof=5Findices (elem, dof=5Findices=5Fp, p=5Fvar); >=20 > const unsigned int n=5Fdofs =3D dof=5Findices.size(); > const unsigned int n=5Fu=5Fdofs =3D dof=5Findices=5Fu.size();=20= > const unsigned int n=5Fv=5Fdofs =3D dof=5Findices=5Fv.size();=20= > const unsigned int n=5Fp=5Fdofs =3D dof=5Findices=5Fp.size(); > =20 > // Compute the elementspecific data for the current > // element. This involves computing the location of the > // quadrature points (q=5Fpoint) and the shape functions > // (phi, dphi) for the current element. > fe=5Fvel>reinit (elem); > fe=5Fpres>reinit (elem); >=20 > // Zero the element matrix and righthand side before > // summing them. We use the resize member here because > // the number of degrees of freedom might have changed from > // the last element. Note that this will be the case if the > // element type is different (i.e. the last element was a > // triangle, now we are on a quadrilateral). > Ke.resize (n=5Fdofs, n=5Fdofs); > Fe.resize (n=5Fdofs); >=20 > // Reposition the submatrices... The idea is this: > // > //     > //  Kuu Kuv Kup   Fu  > // Ke =3D  Kvu Kvv Kvp ; Fe =3D  Fv  > //  Kpu Kpv Kpp   Fp  > //     > // > // The \p DenseSubMatrix.repostition () member takes the > // (row=5Foffset, column=5Foffset, row=5Fsize, column=5Fsize).= > //=20 > // Similarly, the \p DenseSubVector.reposition () member > // takes the (row=5Foffset, row=5Fsize) > Kuu.reposition (u=5Fvar*n=5Fu=5Fdofs, u=5Fvar*n=5Fu=5Fdofs, n=5F= u=5Fdofs, n=5Fu=5Fdofs); > Kuv.reposition (u=5Fvar*n=5Fu=5Fdofs, v=5Fvar*n=5Fu=5Fdofs, n=5F= u=5Fdofs, n=5Fv=5Fdofs); > Kup.reposition (u=5Fvar*n=5Fu=5Fdofs, p=5Fvar*n=5Fu=5Fdofs, n=5F= u=5Fdofs, n=5Fp=5Fdofs); > =20 > Kvu.reposition (v=5Fvar*n=5Fv=5Fdofs, u=5Fvar*n=5Fv=5Fdofs, n=5F= v=5Fdofs, n=5Fu=5Fdofs); > Kvv.reposition (v=5Fvar*n=5Fv=5Fdofs, v=5Fvar*n=5Fv=5Fdofs, n=5F= v=5Fdofs, n=5Fv=5Fdofs); > Kvp.reposition (v=5Fvar*n=5Fv=5Fdofs, p=5Fvar*n=5Fv=5Fdofs, n=5F= v=5Fdofs, n=5Fp=5Fdofs); > =20 > Kpu.reposition (p=5Fvar*n=5Fu=5Fdofs, u=5Fvar*n=5Fu=5Fdofs, n=5F= p=5Fdofs, n=5Fu=5Fdofs); > Kpv.reposition (p=5Fvar*n=5Fu=5Fdofs, v=5Fvar*n=5Fu=5Fdofs, n=5F= p=5Fdofs, n=5Fv=5Fdofs); > Kpp.reposition (p=5Fvar*n=5Fu=5Fdofs, p=5Fvar*n=5Fu=5Fdofs, n=5F= p=5Fdofs, n=5Fp=5Fdofs); >=20 > Fu.reposition (u=5Fvar*n=5Fu=5Fdofs, n=5Fu=5Fdofs); > Fv.reposition (v=5Fvar*n=5Fu=5Fdofs, n=5Fv=5Fdofs); > Fp.reposition (p=5Fvar*n=5Fu=5Fdofs, n=5Fp=5Fdofs); >=20 > // Now we will build the element matrix and righthandside. > // Constructing the RHS requires the solution and its > // gradient from the previous timestep. This must be > // calculated at each quadrature point by summing the > // solution degreeoffreedom values by the appropriate > // weight functions. > for (unsigned int qp=3D0; qp<qrule.n=5Fpoints(); qp++) > =09{ > =09 // Values to hold the solution & its gradient at the previous t= imestep. > =09 Number u =3D 0., u=5Fold =3D 0.; > =09 Number v =3D 0., v=5Fold =3D 0.; > =09 Number p=5Fold =3D 0.; > =09 Gradient grad=5Fu, grad=5Fu=5Fold; > =09 Gradient grad=5Fv, grad=5Fv=5Fold; > =09 =20 > =09 // Compute the velocity & its gradient from the previous timest= ep > =09 // and the old Newton iterate. > =09 for (unsigned int l=3D0; l<n=5Fu=5Fdofs; l++) > =09 { > =09 // From the old timestep: > =09 u=5Fold +=3D phi[l][qp]*navier=5Fstokes=5Fsystem.old=5Fsolu= tion (dof=5Findices=5Fu[l]); > =09 v=5Fold +=3D phi[l][qp]*navier=5Fstokes=5Fsystem.old=5Fsolu= tion (dof=5Findices=5Fv[l]); > =09 grad=5Fu=5Fold.add=5Fscaled (dphi[l][qp],navier=5Fstokes=5F= system.old=5Fsolution (dof=5Findices=5Fu[l])); > =09 grad=5Fv=5Fold.add=5Fscaled (dphi[l][qp],navier=5Fstokes=5F= system.old=5Fsolution (dof=5Findices=5Fv[l])); >=20 > =09 // From the previous Newton iterate: > =09 u +=3D phi[l][qp]*navier=5Fstokes=5Fsystem.current=5Fsoluti= on (dof=5Findices=5Fu[l]);=20 > =09 v +=3D phi[l][qp]*navier=5Fstokes=5Fsystem.current=5Fsoluti= on (dof=5Findices=5Fv[l]); > =09 grad=5Fu.add=5Fscaled (dphi[l][qp],navier=5Fstokes=5Fsystem= .current=5Fsolution (dof=5Findices=5Fu[l])); > =09 grad=5Fv.add=5Fscaled (dphi[l][qp],navier=5Fstokes=5Fsystem= .current=5Fsolution (dof=5Findices=5Fv[l])); > =09 } >=20 > =09 // Compute the old pressure value at this quadrature point. > =09 for (unsigned int l=3D0; l<n=5Fp=5Fdofs; l++) > =09 { > =09 p=5Fold +=3D psi[l][qp]*navier=5Fstokes=5Fsystem.old=5Fsolu= tion (dof=5Findices=5Fp[l]); > =09 } >=20 > =09 // Definitions for convenience. It is sometimes simpler to do = a > =09 // dot product if you have the full vector at your disposal. > =09 const NumberVectorValue U=5Fold (u=5Fold, v=5Fold); > =09 const NumberVectorValue U (u, v); > =09 const Number u=5Fx =3D grad=5Fu(0); > =09 const Number u=5Fy =3D grad=5Fu(1); > =09 const Number v=5Fx =3D grad=5Fv(0); > =09 const Number v=5Fy =3D grad=5Fv(1); > =09 =20 > =09 // First, an iloop over the velocity degrees of freedom. > =09 // We know that n=5Fu=5Fdofs =3D=3D n=5Fv=5Fdofs so we can comp= ute contributions > =09 // for both at the same time. > =09 for (unsigned int i=3D0; i<n=5Fu=5Fdofs; i++) > =09 { > =09 Fu(i) +=3D JxW[qp]*(u=5Fold*phi[i][qp]  = // massmatrix term=20 > =09=09=09=09(1.theta)*dt*(U=5Fold*grad=5Fu=5Fold)*phi[i][qp] + // c= onvection term > =09=09=09=09(1.theta)*dt*p=5Fold*dphi[i][qp](0)  // press= ure term on rhs > =09=09=09=09(1.theta)*dt*(grad=5Fu=5Fold*dphi[i][qp]) + // dif= fusion term on rhs > =09=09=09=09theta*dt*(U*grad=5Fu)*phi[i][qp]); // Newto= n term >=20 > =09=09 > =09 Fv(i) +=3D JxW[qp]*(v=5Fold*phi[i][qp]  = // massmatrix term > =09=09=09=09(1.theta)*dt*(U=5Fold*grad=5Fv=5Fold)*phi[i][qp] + // = convection term > =09=09=09=09(1.theta)*dt*p=5Fold*dphi[i][qp](1)  // pres= sure term on rhs > =09=09=09=09(1.theta)*dt*(grad=5Fv=5Fold*dphi[i][qp]) + // di= ffusion term on rhs > =09=09=09=09theta*dt*(U*grad=5Fv)*phi[i][qp]); // Newt= on term > =09 =20 >=20 > =09 // Note that the Fp block is identically zero unless we are= using > =09 // some kind of artificial compressibility scheme... >=20 > =09 // Matrix contributions for the uu and vv couplings. > =09 for (unsigned int j=3D0; j<n=5Fu=5Fdofs; j++) > =09=09{ > =09=09 Kuu(i,j) +=3D JxW[qp]*(phi[i][qp]*phi[j][qp] + = // mass matrix term > =09=09=09=09 theta*dt*(dphi[i][qp]*dphi[j][qp]) + // diffusi= on term > =09=09=09=09 theta*dt*(U*dphi[j][qp])*phi[i][qp] + // convect= ion term > =09=09=09=09 theta*dt*u=5Fx*phi[i][qp]*phi[j][qp]); // Newto= n term >=20 > =09=09 Kuv(i,j) +=3D JxW[qp]*theta*dt*u=5Fy*phi[i][qp]*phi[j][qp]; = // Newton term > =09=09 =20 > =09=09 Kvv(i,j) +=3D JxW[qp]*(phi[i][qp]*phi[j][qp] + = // mass matrix term > =09=09=09=09 theta*dt*(dphi[i][qp]*dphi[j][qp]) + // diffusi= on term > =09=09=09=09 theta*dt*(U*dphi[j][qp])*phi[i][qp] + // convect= ion term > =09=09=09=09 theta*dt*v=5Fy*phi[i][qp]*phi[j][qp]); // Newto= n term >=20 > =09=09 Kvu(i,j) +=3D JxW[qp]*theta*dt*v=5Fx*phi[i][qp]*phi[j][qp]; = // Newton term > =09=09} >=20 > =09 // Matrix contributions for the up and vp couplings. > =09 for (unsigned int j=3D0; j<n=5Fp=5Fdofs; j++) > =09=09{ > =09=09 Kup(i,j) +=3D JxW[qp]*(theta*dt*psi[j][qp]*dphi[i][qp](0));= > =09=09 Kvp(i,j) +=3D JxW[qp]*(theta*dt*psi[j][qp]*dphi[i][qp](1));= > =09=09} > =09 } >=20 > =09 // Now an iloop over the pressure degrees of freedom. This co= de computes > =09 // the matrix entries due to the continuity equation. Note: To= maintain a > =09 // symmetric matrix, we may (or may not) multiply the continuit= y equation by > =09 // negative one. Here we do not. > =09 for (unsigned int i=3D0; i<n=5Fp=5Fdofs; i++) > =09 { > =09 for (unsigned int j=3D0; j<n=5Fu=5Fdofs; j++) > =09=09{ > =09=09 Kpu(i,j) +=3D JxW[qp]*psi[i][qp]*dphi[j][qp](0); > =09=09 Kpv(i,j) +=3D JxW[qp]*psi[i][qp]*dphi[j][qp](1); > =09=09} > =09 } > =09} // end of the quadrature point qploop >=20 > =20 > // At this point the interior element integration has > // been completed. However, we have not yet addressed > // boundary conditions. For this example we will only > // consider simple Dirichlet boundary conditions imposed > // via the penalty method. The penalty method used here > // is equivalent (for Lagrange basis functions) to lumping > // the matrix resulting from the L2 projection penalty > // approach introduced in example 3. > { > =09// The following loops over the sides of the element. > =09// If the element has no neighbor on a side then that > =09// side MUST live on a boundary of the domain. > =09for (unsigned int s=3D0; s<elem>n=5Fsides(); s++) > =09 if (elem>neighbor(s) =3D=3D NULL) > =09 { > =09 // Get the boundary ID for side 's'. > =09 // These are set internally by build=5Fsquare(). > =09 // 0=3Dbottom > =09 // 1=3Dright > =09 // 2=3Dtop > =09 // 3=3Dleft > =09 short int bc=5Fid =3D mesh.boundary=5Finfo>boundary=5Fid (= elem,s); > =09 if (bc=5Fid=3D=3DBoundaryInfo::invalid=5Fid) > =09=09 error(); >=20 > =09 =20 > =09 AutoPtr<Elem> side (elem>build=5Fside(s)); > =09 =09 =20 > =09 // Loop over the nodes on the side. > =09 for (unsigned int ns=3D0; ns<side>n=5Fnodes(); ns++) > =09=09{ > =09=09 // The penalty value. \f$ \frac{1}{\epsilon \f$ > =09=09 const Real penalty =3D 1.e10; > =09=09 =20 > =09=09 // Get the boundary values. > =09=09 =20 > =09=09 // Set u =3D 1 on the top boundary, 0 everywhere else > =09=09 const Real u=5Fvalue =3D (bc=5Fid=3D=3D2) =3F 1. : 0.; > =09=09 =20 > =09=09 // Set v =3D 0 everywhere > =09=09 const Real v=5Fvalue =3D 0.; > =09=09 =20 > =09=09 // Find the node on the element matching this node on > =09=09 // the side. That defined where in the element matrix > =09=09 // the boundary condition will be applied. > =09=09 for (unsigned int n=3D0; n<elem>n=5Fnodes(); n++) > =09=09 if (elem>node(n) =3D=3D side>node(ns)) > =09=09 { > =09=09=09// Matrix contribution. > =09=09=09Kuu(n,n) +=3D penalty; > =09=09=09Kvv(n,n) +=3D penalty; > =09=09 =09=09 =20 > =09=09=09// Righthandside contribution. > =09=09=09Fu(n) +=3D penalty*u=5Fvalue; > =09=09=09Fv(n) +=3D penalty*v=5Fvalue; > =09=09 } > =09=09} // end face node loop=09 =20 > =09 } // end if (elem>neighbor(side) =3D=3D NULL) >=20 > =09/* qqq begin */ >=20 > =09for (unsigned int c=3D0; c<elem>n=5Fnodes(); c++) > =09 if (elem>node(c) =3D=3D 0) > =09 { > =09 const Real penalty =3D 1.e10; > =09 Kpp(c,c) +=3D penalty; > =09 Fp(c) +=3D penalty*100.0; > =09 } > =09 > =09/* qqq end */ >=20 > } // end boundary condition section=09 =20 > =20 > /* qqq begin */ > if(fabs(elem>centroid()(0)0.5)<=3D0.4 && fabs(elem>centroid= ()(1)0.5)<=3D0.4) > =09{ > =09 const Real penalty =3D 1.e7; > =09 for (unsigned int i=3D0; i<n=5Fu=5Fdofs; i++) > =09 { > =09 Kuu(i,i) +=3D penalty; > =09 Kvv(i,i) +=3D penalty; > =09 } > =09} > /* qqq end */ >=20 > // The element matrix and righthandside are now built > // for this element. Add them to the global matrix and > // righthandside vector. The \p PetscMatrix::add=5Fmatrix()= > // and \p PetscVector::add=5Fvector() members do this for us. > navier=5Fstokes=5Fsystem.matrix>add=5Fmatrix (Ke, dof=5Findic= es); > navier=5Fstokes=5Fsystem.rhs>add=5Fvector (Fe, dof=5Findic= es); > } // end of element loop > =20 > // That's it. > return; > } 
From: <tim@ce...>  20070116 08:14:30

Dear Roy and John, thank you for your nearly simultaneous answer. (Did you use a stopwatch to get that right?) (: Actually, constraining away the onedimensional null space did not cure the problem in my case. Please find attached a modified version of my test program that includes John's code snippet. Again, all the (very little) places where the program differs from ex13 are marked with "qqq". The behaviour is the same as before, i.e. the program works well without obstacle (rough and fine grid), and it works well with obstacle on the rough grid (20x20), but it crashes with obstacle on the fine grid (100x100). I pinned the pressure to 100.0 at node 0, but pinning it to 0.0 instead does not change the problem. It would be nice if you at least tested whether this phenomenon also happens with your installation, in order to exclude version conflicts. Thank you very much in advance, Tim 