From: Rolf Mueller <rolfm@mi...>  20050310 08:14:06

Dear libmesh developers and users, I would like to obtain complex (i.e., real and imaginary) field values for a problem with an unbounded domain. I have been looking into the libmesh example no. 6 (ex6) in order to see how to use libmesh's infinite elements. I have also read the paper on the wave envelope elements by Astley et al. (J. Acoust. Soc. Am. 103, 4963, 1998) and perused the paper by Dreyer et al. (Int. J. Numer. Meth. Engng. 58, 933953, 2003). In the paper by Astley et al., the contributions of the infinite elements to the system matrix appear to be complex; Eq. 33 of the paper states they should be: K+jkCk^2M, where K,C,M are the stiffness, damping, and mass matrix respectively, k is the wavenumber (w/c), and j the imaginary unit. If I look at examples/ex6/ex6.C, I see that the element contributions to the system matrix are assembled as Ke.add(1./speed , Ce); Ke.add(1./(speed*speed), Me); The real part of the damping matrix Ce seems to be  in general  nonzero, whereas the imaginary part seems to be always zero (I have configured libmesh with the "enablecomplex"flag). So it seems that Ce is just added to a real system matrix. In addition, the sign of Me appears to be the opposite of what I would have expected from Eq. 33 in Astley et al. ex6 seems to compute a complex result, the "System "Wave" Solution Vector" written to the file "eqn_sys.dat" is twice as long as the number of nodes. However, the values (real and imaginary) seem to be always zero, even if not all nodes are part of infinite elements (I have been using meshe sizes up to 20x20x20). Could somebody please help me understand what is being done here and how I can get a physically meaningful result for the (finite) nodes?  I am thankful for any input. Kind regards, Rolf 
From: Steffen Petersen <steffen.petersen@tu...>  20050318 00:18:22

Sorry for the late reply, I have been on a short vacation + conference trip. Rolf Mueller wrote: > Dear libmesh developers and users, > > I would like to obtain complex (i.e., real and imaginary) field values > for a problem with an unbounded domain. > > I have been looking into the libmesh example no. 6 (ex6) in order to > see how to use libmesh's infinite elements. I have also read the paper > on the wave envelope elements by Astley et al. > (J. Acoust. Soc. Am. 103, 4963, 1998) and perused the paper by Dreyer > et al. (Int. J. Numer. Meth. Engng. 58, 933953, 2003). > > In the paper by Astley et al., the contributions of the infinite > elements to the system matrix appear to be complex; Eq. 33 of the > paper states they should be: > > K+jkCk^2M, > > where K,C,M are the stiffness, damping, and mass matrix respectively, > k is the wavenumber (w/c), and j the imaginary unit. > > If I look at examples/ex6/ex6.C, I see that the element contributions > to the system matrix are assembled as > > Ke.add(1./speed , Ce); > Ke.add(1./(speed*speed), Me); > The real part of the damping matrix Ce seems to be  in general  > nonzero, whereas the imaginary part seems to be always zero (I have > configured libmesh with the "enablecomplex"flag). So it seems that > Ce is just added to a real system matrix. In addition, the sign of Me > appears to be the opposite of what I would have expected from Eq. 33 > in Astley et al. The reason for the confusions is that in example 6 the matrix assembly is done according to transient simulations (cf. part II of the Astley paper; however, no time stepping scheme is accomplished here). To obtain the system matrix for solving the Helmholtz equation (as is done in the papers you mention) you may simply change the factors for K, C and M, e.g.: const Complex scale_damping ( 0., k ); const Complex scale_mass (k*k, 0.); Ke.add(scale_damping, Ce); Ke.add(scale_mass , Me); If you have compiled the code with enablecomplex, the system matrices, rhs and solution will always be complex (also in the example programs covering real valued problems). > ex6 seems to compute a complex result, the "System "Wave" Solution > Vector" written to the file "eqn_sys.dat" is twice as long as the > number of nodes. However, the values (real and imaginary) seem to be > always zero, even if not all nodes are part of infinite elements (I > have been using meshe sizes up to 20x20x20). I just ran the example program and could not reproduce your problems. The vector contains complex data due to the configuration, where the imaginary part is zero since matrix and rhs are real. Here are the first lines of the eqn_sys.dat file I get: 1 # No. of Equation Systems Wave # Name, System No. 0 LinearImplicit # Type, System No. 0 1 # No. of Variables in System "Wave" p # Name, Variable No. 0, System "Wave" 1 # Approximation Order, Variable "p", System "Wave" 3 # Radial Approximation Order, Variable "p", System "Wave" 0 # FE Family, Variable "p", System "Wave" 12 # Radial FE Family, Variable "p", System "Wave" 0 # Infinite Mapping Type, Variable "p", System "Wave" 0 # No. of Additional Vectors, System "Wave" 419 # vector length x 2 (complex) 8.491741e03 0.000000e+00 1.204387e02 0.000000e+00 1.606009e02 0.000000e+00 1.204387e02 0.000000e+00 8.491741e03 0.000000e+00 1.204387e02 0.000000e+00 2.047754e02 0.000000e+00 1.793855e02 0. > Could somebody please help me understand what is being done here and > how I can get a physically meaningful result for the (finite) nodes?  The nodal data within the finite element region should be physically meaningful. Only for evaluation of results in the infinte domain special treatment is necessary. Steffen > I am thankful for any input. > > Kind regards, > Rolf > > >  > 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: Rolf Mueller <rolfm@mi...>  20050323 11:03:11

Dear Steffen, thank you so much for your kind help. I appreciate it a lot. On Fri, Mar 18, 2005 at 01:18:00AM +0100, Steffen Petersen wrote: > The reason for the confusions is that in example 6 the matrix assembly > is done according to transient simulations > (cf. part II of the Astley paper; however, no time stepping scheme is > accomplished here). > To obtain the system matrix for solving the Helmholtz equation (as is > done in the papers you mention) you may simply change the factors for > K, C and M, e.g.: > > const Complex scale_damping ( 0., k ); > const Complex scale_mass (k*k, 0.); > > Ke.add(scale_damping, Ce); > Ke.add(scale_mass , Me); Excellent!  Many thanks again for helping me out of my confusion here. > I just ran the example program and could not reproduce your > problems. The vector contains complex data due to the > configuration, where the imaginary part is zero since matrix and rhs > are real. Here are the first lines of the eqn_sys.dat file I get: I have now compiled a fresh version of libmesh0.4.3_rc2 with the configuration options "enablelaspack disablempi disablepetsc". This gives me the same numbers you are reporting for the original timedomain system matrix. If I set up the system matrix for the frequencydomain Helmholtz equation in the way you suggest, I also get both nonzero real and imaginary values. > The nodal data within the finite element region should be physically > meaningful. Only for evaluation of results in the infinte domain > special treatment is necessary. That's good to know. However, I am still not sure if I interpret the results correctly: In ex6, "a unit load" is applied "at the node located at (0,0,0)". This seems to be done by adding 1 to the rhs vector at the position of this node. I have looked at the results as written out by "GMVIO(mesh).write_equation_systems". Shouldn't I see something like 1/(4*pi*r)*exp(j*k*r), were r=sqrt(x.^2+y.^2+z.^2) and k is the wavenumber? I have chosen k=1 and 10 or 20 elements along each dimension in the cube [1 1][1 1][1 1]. What I get are positive and negative numbers in the order of magnitude of 10^42 to 10^45 or so for both the real and imaginary parts when using 10 elements along each dimension. For 20 elements, the magnitude of the numbers appears to be much smaller (10^31  10^36). Besides the values, the geometry of the solution is also different from what I would have expected (spherical symmetry). There is only a smear from near the center into the direction of one cube corner.  Do you happen to see where the confusion now lies? Thank you so much once again. Rolf  Rolf Mueller http://www.mip.sdu.dk/~rolfm Maersk Institute, University of Southern Denmark Campusvej 55, DK5230 Odense M, Denmark Phone ++45 6550 3655 Fax. ++45 6615 7697 
From: Steffen Petersen <steffen.petersen@tu...>  20050324 17:47:01

> That's good to know. However, I am still not sure if I > interpret the results correctly: In ex6, "a unit load" is > applied "at the node located at (0,0,0)". This seems to be > done by adding 1 to the rhs vector at the position of this > node. I have looked at the results as written out by > "GMVIO(mesh).write_equation_systems". Shouldn't I see > something like 1/(4*pi*r)*exp(j*k*r), were > r=sqrt(x.^2+y.^2+z.^2) and k is the wavenumber? I have chosen > k=1 and 10 or 20 elements along each dimension in the cube > [1 1][1 1][1 1]. What I get are positive and negative > numbers in the order of magnitude of 10^42 to 10^45 or so for > both the real and imaginary parts when using 10 elements > along each dimension. For 20 elements, the magnitude of the > numbers appears to be much smaller (10^31  10^36). Besides > the values, the geometry of the solution is also different > from what I would have expected (spherical symmetry). There > is only a smear from near the center into the direction of > one cube corner.  Do you happen to see where the confusion now lies? The problem is somehow related to laspack. When using PETSc everything works fine and I get reasonable results. With disablepetsc I also obtain the (definitely) wrong results you report. I have to admit that I have not used laspack with complex numbers for quite some time. Hm, I will see if I can find the problem with the interface. Steffen > Thank you so much once again. > Rolf > >  > Rolf Mueller http://www.mip.sdu.dk/~rolfm > Maersk Institute, University of Southern Denmark > Campusvej 55, DK5230 Odense M, Denmark > Phone ++45 6550 3655 Fax. ++45 6615 7697 > > >  > This SF.net email is sponsored by: 2005 Windows Mobile > Application Contest Submit applications for Windows > Mobile(tm)based Pocket PCs or Smartphones for the chance to > win $25,000 and application distribution. Enter today at > http://ads.osdn.com/?ad_id=6882&alloc_id=15148&op=click > _______________________________________________ > Libmeshusers mailing list > Libmeshusers@... > https://lists.sourceforge.net/lists/listinfo/libmeshusers > 
From: Rolf Mueller <rolfm@mi...>  20050326 20:16:54

Dear Steffen, thank you so much again for your kind help. On Thu, Mar 24, 2005 at 06:46:39PM +0100, Steffen Petersen wrote: > The problem is somehow related to laspack. When using PETSc > everything works fine and I get reasonable results. OK, what an irony: I had chosen laspack, because I had  wrongly  assumed that it would be the least errorprone. Now I am using PETSc (configuration "enablepetsc enablecomplex enablempi") and I think it works: The spherical symmetry of the solution looks very good in general and was perfect where I tested it (along the xaxis). Here are some example results for the real part of the solution along the xaxis: x ex6 cos(x)/(4*pi*x) 0 6.5665 Inf 0.0667 0.6505 1.1910 0.1333 0.6000 0.5915 0.2000 0.3743 0.3900 0.2667 0.2840 0.2879 0.3333 0.2231 0.2256 0.4000 0.1818 0.1832 0.4667 0.1514 0.1523 0.5333 0.1278 0.1285 0.6000 0.1090 0.1095 0.6667 0.0934 0.0938 0.7333 0.0803 0.0806 0.8000 0.0691 0.0693 0.8667 0.0592 0.0594 0.9333 0.0506 0.0507 1.0000 0.0214 0.0430 So it seems that  away from the singularity  the accuracy of the solution is pretty good, but it deteriorates at the node which is shared with the infinite element (at x=1).  Is that a necessity, BTW? Thank you once more. Happy Easter, Rolf 
From: Steffen Petersen <steffen.petersen@tu...>  20050329 10:02:43

> > OK, what an irony: I had chosen laspack, because I had  > wrongly  assumed that it would be the least errorprone. Now > I am using PETSc (configuration "enablepetsc > enablecomplex enablempi") and I think it works: The > spherical symmetry of the solution looks very good in general > and was perfect where I tested it (along the xaxis). > > Here are some example results for the real part of the > solution along the xaxis: > > x ex6 cos(x)/(4*pi*x) > 0 6.5665 Inf > 0.0667 0.6505 1.1910 > 0.1333 0.6000 0.5915 > 0.2000 0.3743 0.3900 > 0.2667 0.2840 0.2879 > 0.3333 0.2231 0.2256 > 0.4000 0.1818 0.1832 > 0.4667 0.1514 0.1523 > 0.5333 0.1278 0.1285 > 0.6000 0.1090 0.1095 > 0.6667 0.0934 0.0938 > 0.7333 0.0803 0.0806 > 0.8000 0.0691 0.0693 > 0.8667 0.0592 0.0594 > 0.9333 0.0506 0.0507 > 1.0000 0.0214 0.0430 > > So it seems that  away from the singularity  the accuracy > of the solution is pretty good, but it deteriorates at the > node which is shared with the infinite element (at x=1).  Is > that a necessity, BTW? The main reason for these inaccuracies is probably the box shaped infinite element interface (envelope). A spherical envelope should lead to more accurate results. Note that you can also increase the radial order of the infinite elements to increase the accuracy of your computations. Regarding example 6, this can be done changing line 163 FEType fe_type(FIRST); to FEType fe_type(FIRST, LAGRANGE, requested_ifem_order); The default radial order is THIRD. The infinite elements are implemented up to order EIGHTTEENTH, which is perhaps solely of academic use. Steffen 
From: Mark Blome <mblome@au...>  20050329 13:47:18

Hi Rolf, > > > > So it seems that  away from the singularity  the accuracy > > of the solution is pretty good, but it deteriorates at the > > node which is shared with the infinite element (at x=3D1).  Is > > that a necessity, BTW? > I noticed the same in my calculations. It seems like the solution becomes=20 inaccurate close to the nodes that are shared with the infinite elements bu= t=20 =2D astonishingly becomes accurate again at some distance to these nodes. = I=20 tried different settings for the infinite elements radial order (and=20 different polynomials), but that doesnt improve things. Also I tried, as=20 Steffen suggested, a spherical mesh (please see the two pictures attached,= =20 they show scalar cut planes through 3D calculations, the border of the mesh= =20 is where the distortions occur, in this case there are no infinite elements= =20 attached at the upper surface ), but the problem remains.=20 I'll let you know if I find out more about whats going wrong, Mark Am Dienstag, 29. M=E4rz 2005 12.02 schrieb Steffen Petersen: > > OK, what an irony: I had chosen laspack, because I had  > > wrongly  assumed that it would be the least errorprone. Now > > I am using PETSc (configuration "enablepetsc > > enablecomplex enablempi") and I think it works: The > > spherical symmetry of the solution looks very good in general > > and was perfect where I tested it (along the xaxis). > > > > Here are some example results for the real part of the > > solution along the xaxis: > > > > x ex6 cos(x)/(4*pi*x) > > 0 6.5665 Inf > > 0.0667 0.6505 1.1910 > > 0.1333 0.6000 0.5915 > > 0.2000 0.3743 0.3900 > > 0.2667 0.2840 0.2879 > > 0.3333 0.2231 0.2256 > > 0.4000 0.1818 0.1832 > > 0.4667 0.1514 0.1523 > > 0.5333 0.1278 0.1285 > > 0.6000 0.1090 0.1095 > > 0.6667 0.0934 0.0938 > > 0.7333 0.0803 0.0806 > > 0.8000 0.0691 0.0693 > > 0.8667 0.0592 0.0594 > > 0.9333 0.0506 0.0507 > > 1.0000 0.0214 0.0430 > > > > So it seems that  away from the singularity  the accuracy > > of the solution is pretty good, but it deteriorates at the > > node which is shared with the infinite element (at x=3D1).  Is > > that a necessity, BTW? > > The main reason for these inaccuracies is probably the > box shaped infinite element interface (envelope). > A spherical envelope should lead to more accurate results. > Note that you can also increase the radial order of the > infinite elements to increase the accuracy of your computations. > Regarding example 6, this can be done changing line 163 > FEType fe_type(FIRST); to FEType fe_type(FIRST, LAGRANGE, > requested_ifem_order); > The default radial order is THIRD. The infinite elements are > implemented up to order EIGHTTEENTH, which is perhaps > solely of academic use. > > Steffen > > > >  > 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: Steffen Petersen <steffen.petersen@tu...>  20050401 12:47:08

Hi Rolf, hi Mark, I have run example 6 with the changes Rolf has made for the matrix assembly here are the nodal results along the xaxis=20 I get using 10 elements in each direction: x Re(p) 0.200000 0.209618 0.400000 0.185383 0.600000 0.103573 0.800000 0.0673304 1.000000 0.0415307 Away from the singularity, this seems to be in agreement with cos(x)/(4*pi*x) (considering the unfavorable shape of the envelope). I have also evaluated results on a fieldpoint line that crosses the envelope (see results below). At this, I also get sufficiently accurat results close to the envelope (x=3D1). Hence, the problem may stem from the evaluation of results. How do you proceed to obtain your results? Steffen x Re(p) 0.2 0.209618=09 0.266667 0.20154=09 0.333333 0.193461=09 0.4 0.185383=09 0.466667 0.158113=09 0.533333 0.130843=09 0.6 0.103573=09 0.666667 0.0914922=09 0.733333 0.0794113=09 0.8 0.0673304=09 0.866667 0.0587305=09 0.933333 0.0501306=09 1 0.0415307=09 1.06667 0.0348121=09 1.13333 0.0287062=09 1.2 0.0231339=09 1.26667 0.0180338=09 1.33333 0.0133562=09 1.4 0.00906065=09 1.46667 0.00511327=09 1.53333 0.00148564=09 1.6 0.00184646=09 1.66667 0.00490379=09 1.73333 0.00770433=09 1.8 0.0102637=09 1.86667 0.0125958=09 1.93333 0.0147126=09 2 0.0166252=09 Steffen >=20 > Hi Rolf, > > > > > > So it seems that  away from the singularity  the=20 > accuracy of the=20 > > > solution is pretty good, but it deteriorates at the node which is=20 > > > shared with the infinite element (at x=3D1).  Is that a = necessity,=20 > > > BTW? > > > I noticed the same in my calculations. It seems like the=20 > solution becomes=20 > inaccurate close to the nodes that are shared with the=20 > infinite elements but=20 >  astonishingly becomes accurate again at some distance to=20 > these nodes. I=20 > tried different settings for the infinite elements radial order (and=20 > different polynomials), but that doesnt improve things. Also=20 > I tried, as=20 > Steffen suggested, a spherical mesh (please see the two=20 > pictures attached,=20 > they show scalar cut planes through 3D calculations, the=20 > border of the mesh=20 > is where the distortions occur, in this case there are no=20 > infinite elements=20 > attached at the upper surface ), but the problem remains.=20 > I'll let you know if I find out more about whats going wrong, >=20 > Mark >=20 >=20 > Am Dienstag, 29. M=E4rz 2005 12.02 schrieb Steffen Petersen: > > > OK, what an irony: I had chosen laspack, because I had =20 > wrongly =20 > > > assumed that it would be the least errorprone. Now I am=20 > using PETSc=20 > > > (configuration "enablepetsc enablecomplex=20 > enablempi") and I=20 > > > think it works: The spherical symmetry of the solution looks very=20 > > > good in general and was perfect where I tested it (along the=20 > > > xaxis). > > > > > > Here are some example results for the real part of the solution=20 > > > along the xaxis: > > > > > > x ex6 cos(x)/(4*pi*x) > > > 0 6.5665 Inf > > > 0.0667 0.6505 1.1910 > > > 0.1333 0.6000 0.5915 > > > 0.2000 0.3743 0.3900 > > > 0.2667 0.2840 0.2879 > > > 0.3333 0.2231 0.2256 > > > 0.4000 0.1818 0.1832 > > > 0.4667 0.1514 0.1523 > > > 0.5333 0.1278 0.1285 > > > 0.6000 0.1090 0.1095 > > > 0.6667 0.0934 0.0938 > > > 0.7333 0.0803 0.0806 > > > 0.8000 0.0691 0.0693 > > > 0.8667 0.0592 0.0594 > > > 0.9333 0.0506 0.0507 > > > 1.0000 0.0214 0.0430 > > > > > > So it seems that  away from the singularity  the=20 > accuracy of the=20 > > > solution is pretty good, but it deteriorates at the node which is=20 > > > shared with the infinite element (at x=3D1).  Is that a = necessity,=20 > > > BTW? > > > > The main reason for these inaccuracies is probably the > > box shaped infinite element interface (envelope). > > A spherical envelope should lead to more accurate results.=20 > Note that=20 > > you can also increase the radial order of the infinite elements to=20 > > increase the accuracy of your computations. Regarding=20 > example 6, this=20 > > can be done changing line 163 FEType fe_type(FIRST); to FEType=20 > > fe_type(FIRST, LAGRANGE, > > requested_ifem_order); > > The default radial order is THIRD. The infinite elements are=20 > > implemented up to order EIGHTTEENTH, which is perhaps solely of=20 > > academic use. > > > > Steffen > > > > > > > >  > > SF email is sponsored by  The IT Product Guide > > Read honest & candid reviews on hundreds of IT Products from real=20 > > users. Discover which products truly live up to the hype. Start=20 > > 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 >=20 