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}
(25) 
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...>  20070115 21:26:54

On Mon, 15 Jan 2007, Travis Austin wrote: > It is probably best to completely replace the init() call with an > init(PetscMatrix <T>* matrix) call. Any thoughts on the cleanest > way to fix this up. Well, we'd need that to be init(SparseMatrix<T>* matrix) instead for consistency with other solver implementations... but then what happens when we don't have a matrix built when we call init? In the PETSc methods, pc() and ksp() both call init(); what do we have them hand in as a matrix? I'd certainly like algebraic multigrid working, but I don't want to inadvertently break anything else in the process. I don't know enough about the PETSc API, but is there really no way to move the offending functions from init() to solve() to avoid such an odd matrix dependency?  Roy 
From: Travis Austin <austint73@ya...>  20070115 20:35:20

Roy and others,=0A=0AI've successfully been able to use HYPRE/BoomerAMG wit= h libMesh. Here is what I did. =0A=0AI compiled libMesh with petsc2.3.2p= 7 and hypre1.1.19b. My first attempt at using =0AHYPRE/BoomerAMG ("pc_ty= pe hypre pc_hypre_type boomeramg") ran into the same hurdle =0Aas found by= Roy. Also, using HYPRE/Euclid ("pc_hypre_type euclid") worked fine.=0A= =0AThe problem was that there is a slightly different initialization proces= s for BoomerAMG within =0AHYPRE and the current set of calls in petsc_linea= r_solver.C works for HYPRE/Euclid but not for HYPRE/BoomerAMG. In particul= ar, it was the init(...) and solve(...) calls that needed minor reworking. = =0AWhat I've done so far is created a new init(...) that takes the matrix = (a PetscMatrix) as an argument =0Aso that I can call KSPSetOperators in ini= t(...). =0A=0AIt is presently called in solve(...) which is too late for H= YPRE/BoomerAMG. HYPRE/BoomerAMG =0Aneeds KSPSetOperators(...) called befor= e KSPSetFromOptions(...). This is all done within init(...) in petsc_linea= r_solver.C. I then commented out KSPSetOperators(...) within solve(...) si= nce it was =0Aalready called.=0A=0AIt is probably best to completely replac= e the init() call with an init(PetscMatrix <T>* matrix) call. =0AAny thoug= hts on the cleanest way to fix this up.=0A=0ACheers,=0ATravis=0A=0A Or= iginal Message =0AFrom: Travis Austin <austint73@...>=0ATo: Roy S= togner <roystgnr@...>; Travis Austin <austint73@...>=0ACc= : libMesh maillist <Libmeshusers@...>=0ASent: Friday, De= cember 15, 2006 8:49:03 AM=0ASubject: Re: [Libmeshusers] Multilevel and mu= ltigrid on libMesh=0A=0ARoy,=0A=0AThat would figure. It looks like I still= need to get hypre installed with my version of petsc. As this is=0Asometh= ing that I want to do for my research (get AMG working) I'll probably be wo= rking on it in the =0Anext several weeks. If anyone else gets it going it = would be nice to hear if there are known hurdles=0Ato get over.=0A=0AAnyway= I assumed that BoomerAMG came installed with hypre by default. I'll check= later.=0A=0ATravis=0A=0A=0A=0A Original Message =0AFrom: Roy Stog= ner <roystgnr@...>=0ATo: Travis Austin <austint73@...>=0A= Cc: libMesh maillist <Libmeshusers@...>=0ASent: Friday, = December 15, 2006 8:13:17 AM=0ASubject: Re: [Libmeshusers] Multilevel and = multigrid on libMesh=0A=0AOn Thu, 14 Dec 2006, Travis Austin wrote:=0A=0A> = Sounds like a great idea and is the direction that we (at the=0A> Universit= y of Auckland) want to take libMesh. But I was figuring on=0A> using Algeb= raic Multigrid (AMG) for all of our multilevel computing=0A> needs. Boomer= AMG from LLNL is available within Petsc (may need to=0A> be downloaded sepa= rately) and easy to use. Is there any particular=0A> reason that you are *= not* using AMG other than better efficiency? I=0A> just wanted to make sur= e you were aware of using AMG as an optimal=0A> solver.=0A=0AI actually did= n't realize that was available; unfortunately even now=0Ait doesn't seem to= be working on my PETSc installation. =0A"pc_type hypre pc_hypre_type euc= lid" works fine, and boomeramg is in=0Athe "help" list of arguments for "= pc_hypre_type", but "pc_type=0Ahypre pc_hypre_type boomeramg" gives me a = segfault. Perhaps I don't=0Ahave the BoomerAMG code installed, and HYPRE o= r PETSc just isn't being=0Asmart about giving me an error message to that e= ffect.=0A=0ARoy=0A=0A= =0ATake Surveys. Earn Cash. Influence the Future of I= T=0AJoin SourceForge.net's Techsay panel and you'll get the chance to share= your=0Aopinions on IT & business topics through brief surveys  and earn c= ash=0Ahttp://www.techsay.com/default.php?page=3Djoin.php&p=3Dsourceforge&CI= D=3DDEVDEV=0A_______________________________________________=0ALibmeshuser= s mailing list=0ALibmeshusers@...=0Ahttps://lists.source= forge.net/lists/listinfo/libmeshusers=0A=0A=0A=0A=0A=0A =0A_______________= _____________________________________________________________________=0ADo = you Yahoo!?=0AEveryone is raving about the allnew Yahoo! Mail beta.=0Ahttp= ://new.mail.yahoo.com=0A=0A= =0ATake Surveys. Earn Cash. Influence the Future o= f IT=0AJoin SourceForge.net's Techsay panel and you'll get the chance to sh= are your=0Aopinions on IT & business topics through brief surveys  and ear= n cash=0Ahttp://www.techsay.com/default.php?page=3Djoin.php&p=3Dsourceforge= &CID=3DDEVDEV=0A_______________________________________________=0ALibmeshu= sers mailing list=0ALibmeshusers@...=0Ahttps://lists.sou= rceforge.net/lists/listinfo/libmeshusers=0A=0A=0A=0A=0A=0A =0A____________= ________________________________________________________________________=0A= Finding fabulous fares is fun. =0ALet Yahoo! FareChase search your favorit= e travel sites to find flight and hotel bargains.=0Ahttp://farechase.yahoo.= com/promogeneric14795097 
From: John Peterson <peterson@cf...>  20070115 17:11:34

./configure ? Seriously though, I don't think LibMesh compiles under windows yet, search through the mailing list, there was one or two postings regarding this in the last few weeks. There are some system functions that libmesh uses which aren't supported under windows I believe. John li pan writes: > Dear developers, > I used to work under linux. Now I'm trying to move the > code into windows. I installed cygwin, and further > installed mpi and petsc. Now I tried to "make" libmesh > in cygwin and got this error: > > $ make > Compiling C++ (in optimized mode) > src/base/dof_map.C... > g++: D:Microsoft: No such file or directory > g++: Visual: No such file or directory > g++: StudioVC98Include: No such file or directory > g++: no input files > /usr/bin/sh: D:QTdevinclude: command not found > /usr/bin/sh: I/tmp/libmesh0.5.0/include/base: No > such file or directory > make: *** [src/base/dof_map.i686pccygwin.opt.o] > Error 127 > > How can I deal with this problem? > > best > > pan 
From: li pan <li76pan@ya...>  20070115 17:06:43

Dear developers, I used to work under linux. Now I'm trying to move the code into windows. I installed cygwin, and further installed mpi and petsc. Now I tried to "make" libmesh in cygwin and got this error: $ make Compiling C++ (in optimized mode) src/base/dof_map.C... g++: D:Microsoft: No such file or directory g++: Visual: No such file or directory g++: StudioVC98Include: No such file or directory g++: no input files /usr/bin/sh: D:QTdevinclude: command not found /usr/bin/sh: I/tmp/libmesh0.5.0/include/base: No such file or directory make: *** [src/base/dof_map.i686pccygwin.opt.o] Error 127 How can I deal with this problem? best pan ____________________________________________________________________________________ Any questions? Get answers on any topic at http://www.Answers.yahoo.com. Try it now. 
From: John Peterson <peterson@cf...>  20070115 15:43:04

One thing that's missing from ex13 is pinning a single value of the pressure. This is a common (at least amongst LibMesh developers) hack for rectifying the nontrivial null space of constant pressures which is associated with the NS equations. Here's a code snippet that can go at the end of your bc routine to fix a single value of the pressure at node "0" (assuming that node 0 is on the boundary.) I will try to get this added to ex13 as well...= John // Pin the pressure to zero const unsigned int pressure=5Fnode=3D0; const Real p=5Fvalue =3D 0.0; for (unsigned int c=3D0; c<elem>n=5Fnodes(); c++) if (elem>node(c) =3D=3D pressure=5Fnode) { =09// here(); =09// std::cout << elem>point(c) << std::endl; =09Kpp(c,c) +=3D penalty; =09Fp(c) +=3D penalty*p=5Fvalue; } Tim Kr=F6ger writes: > Dear libMesh developer team, >=20 > Currently, I'm working on some NavierStokes application with solid=20= > obstacles in the computational domain. My implementation of the=20 > obstacles consists of fixing the velocity to zero on all dofs in the= =20 > respective cells. Please find attached a test program, which is in=20= > fact just a slight modification of ex13; the changes are marked with= =20 > "qqq". (My real application is more complicated, but this suffices = to=20 > demsonstrate my problem.) >=20 > The problem: When using a high resolution grid, no solution is=20 > obtained. Either PETSc crashes with a "Detected zero pivot in LU=20= > factorization" message, or the solver just does not converge. >=20 > Everything works well on a lowresolution grid. Also, everything=20= > works well on a highresolution grid if the obstacle is disabled. >=20 > Of course, another way of implementing the obstacle would be to remo= ve=20 > the obstacle cells from the computational domain, but unfortunately,= =20 > in my real application, I have to solve an additional PDE on the who= le=20 > domain, including the obstacle. Also, the shape and size of the=20 > obstacle vary in time. >=20 > Do you have any ideas what is going wrong and how I can do it better= =3F >=20 > Best Regards, >=20 > Tim >=20 > =20 > Dr. Tim Kroeger > CeVis  Center of Complex Systems and Visualization > University of Bremen > Universitaetsallee 29 (Office 3.13) tim@... > D28359 Bremen Phone +494212187710 > Germany Fax +494212184236/* = $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) > } // 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; > } > =  > Take Surveys. Earn Cash. Influence the Future of IT > Join SourceForge.net's Techsay panel and you'll get the chance to sh= are your > opinions on IT & business topics through brief surveys  and earn ca= sh > http://www.techsay.com/default.php=3Fpage=3Djoin.php&p=3Dsourceforge= &CID=3DDEVDEV=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F= =5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F= =5F=5F=5F > Libmeshusers mailing list > Libmeshusers@... > https://lists.sourceforge.net/lists/listinfo/libmeshusers 
From: Roy Stogner <roystgnr@ic...>  20070115 15:42:35

On Mon, 15 Jan 2007, Tim Kr=F6ger wrote: > Do you have any ideas what is going wrong and how I can do it better? This might just be the second question in a week whose answer is "we don't constrain away the constant pressure solution in ex13". Unless you've pinned a free pressure node along with all your obstacle velocity nodes, your resulting matrix will be one rank deficient, and it's certainly possible that other changes you make to the problem will render the solver (or much more likely, preconditioner) unable to handle that. You might try playing around with different preconditioners before modifying your system, though.  Roy Stogner 
From: <tim@ce...>  20070115 13:43:19

Dear libMesh developer team, Currently, I'm working on some NavierStokes application with solid obstacles in the computational domain. My implementation of the obstacles consists of fixing the velocity to zero on all dofs in the respective cells. Please find attached a test program, which is in fact just a slight modification of ex13; the changes are marked with "qqq". (My real application is more complicated, but this suffices to demsonstrate my problem.) The problem: When using a high resolution grid, no solution is obtained. Either PETSc crashes with a "Detected zero pivot in LU factorization" message, or the solver just does not converge. Everything works well on a lowresolution grid. Also, everything works well on a highresolution grid if the obstacle is disabled. Of course, another way of implementing the obstacle would be to remove the obstacle cells from the computational domain, but unfortunately, in my real application, I have to solve an additional PDE on the whole domain, including the obstacle. Also, the shape and size of the obstacle vary in time. Do you have any ideas what is going wrong and how I can do it better? Best Regards, Tim  Dr. Tim Kroeger CeVis  Center of Complex Systems and Visualization University of Bremen Universitaetsallee 29 (Office 3.13) tim@... D28359 Bremen Phone +494212187710 Germany Fax +494212184236 