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}

_{Jun}

_{Jul}

_{Aug}

_{Sep}

_{Oct}

_{Nov}

_{Dec}

From: Paul T. Bauman <ptbauman@gm...>  20150430 23:49:35

On Thu, Apr 30, 2015 at 7:43 PM, Hafez Asgharzadeh <hafezasg@...> wrote: > > I am wondering how i can get each elements nodes location (X,Y,Z) http://libmesh.github.io/doxygen/classlibMesh_1_1TypeVector.html#a608985f52c18676847e7acb6623beb24 > neighbor elements http://libmesh.github.io/doxygen/classlibMesh_1_1Elem.html#a06576c749e62f4a5962e112e5441c31b > global indices? > http://libmesh.github.io/doxygen/classlibMesh_1_1Elem.html#a7cffa0aae912cc004972c9ea34a22ad4 > Also I am new at Libmesh, Here is the link to the documentation: http://libmesh.github.io/doxygen/classes.html You will find it very useful and you will find folks more willing to help if it's clear you've tried to examine the documentation. > i am wondering how i can use gdb debugger on > libmesh? > Just like any other program. E.g., http://www.unknownroad.com/rtfm/gdbtut/ 
From: Hafez Asgharzadeh <hafezasg@bu...>  20150430 23:44:03

Hi, I am working with 3nodes element. I am wondering how i can get each elements nodes location (X,Y,Z) and also neighbor elements global indices? Also I am new at Libmesh, i am wondering how i can use gdb debugger on libmesh? In advance i appreciate your kind attention, Bests, Hafez 
From: Hafez Asgharzadeh <hafezasg@bu...>  20150429 21:56:43

Hi, I am working on *2D LaplaceYoung Problem Using Nonlinear Solvers.* I am wondering how i can use matrix free snes method instead of computing jacobian? Bests, Hafez 
From: Xujun Zhao <xzhao99@gm...>  20150429 01:34:13

Thank you, David. I will try. Xujun On Tue, Apr 28, 2015 at 5:21 PM, David Knezevic <david.knezevic@...> wrote: > On Tue, Apr 28, 2015 at 6:09 PM, Xujun Zhao <xzhao99@...> wrote: > >> Hi folks, >> >> In my problem, I would like to assemble my global matrix and rhs vector >> separately, because at each time step, my global matrix doesn't change. >> What I need to do is to simply update my rhs vector and resolve it with my >> solver. Is there a way to do this in libMesh? I noticed when we >> attach_assemble_function, both the matrix and rhs vector are >> reconstructed. >> > > > Have a look at System::assemble_before_solve. If you set that to false > then you can call your own assembly functions whenever you like, and the > assembly function won't automatically get called when you call solve. This > way you could implement an assembly function that has bool arguments > assemble_matrix and assemble_rhs and call this whenever you need to with > the bools set appropriately. > > David > > > 
From: Hafez Asgharzadeh <hafezasg@bu...>  20150429 00:10:23

Hi, I am willing to do triangular thin shell element by using TransientNonlinearImplicit System. Assume that this is my equation: Fint(u)+M u"+C u'=Fext where Fint(u)=K(u) S(u), K is a matrix and S is the stress vector. ' is a time derivatives, which i want to discritize it by Newmark method. And for solving this nonlinear equation, i want to use martrxifree method from SNES in petsc. Is there any similar code or example for doing this in Libmesh? which classes can i use from libmesh to ease this procedure? Thanks for your kind attention in advance. Bests, Hafez 
From: Hafez Asgharzadeh <hafezasg@bu...>  20150429 00:06:27

Hi, I am willing to do triangular thin shell element by using TransientNonlinearImplicit System. Assume that this is my equation: Fint(u)+M u"+C u'=Fext where Fint(u)=K(u) S(u), K is a matrix and S is the stress vector. ' is a time derivatives, which i want to discritize it by Newmark method. And for solving this nonlinear equation, i want to use martrxifree method from SNES in petsc. Is there any similar code or example for doing this in Libmesh? which classes can i use from libmesh to ease this procedure? Thanks for your kind attention in advance. Bests, Hafez 
From: David Knezevic <david.knezevic@ak...>  20150428 22:22:02

On Tue, Apr 28, 2015 at 6:09 PM, Xujun Zhao <xzhao99@...> wrote: > Hi folks, > > In my problem, I would like to assemble my global matrix and rhs vector > separately, because at each time step, my global matrix doesn't change. > What I need to do is to simply update my rhs vector and resolve it with my > solver. Is there a way to do this in libMesh? I noticed when we > attach_assemble_function, both the matrix and rhs vector are reconstructed. > Have a look at System::assemble_before_solve. If you set that to false then you can call your own assembly functions whenever you like, and the assembly function won't automatically get called when you call solve. This way you could implement an assembly function that has bool arguments assemble_matrix and assemble_rhs and call this whenever you need to with the bools set appropriately. David 
From: Xujun Zhao <xzhao99@gm...>  20150428 22:09:07

Hi folks, In my problem, I would like to assemble my global matrix and rhs vector separately, because at each time step, my global matrix doesn't change. What I need to do is to simply update my rhs vector and resolve it with my solver. Is there a way to do this in libMesh? I noticed when we attach_assemble_function, both the matrix and rhs vector are reconstructed. Thank you very much. Xujun 
From: Zhang <zyzhang@nu...>  20150428 02:24:43

Dear Vikram, Thank you for the link. But I think there is mistakes in the first relation there. When we write the first equation, generally we consider the Newton's second law applied to the fluid volume, not the solid. Since conventionally \vec{n} points outwards this control volume, then actually the pressure force applied on the surface of fluid control volume by possible solid is minus P other than a positve one ,which then points inside the volume. I think this is the basis of most CFD packages. So in that link most sign of P terms should be changed except the last one, because that is the pressure force exerted by fluid through the surface. Zhenyu 原始邮件 发件人: "Vikram Garg" <vikram.v.garg@...> 发送时间: 20150428 05:39:05 (星期二) 收件人: Zhang <zyzhang@...> 抄送: libmeshusers <libmeshusers@...> 主题: Re: [Libmeshusers] calculation of aerodynamic forces and moments Hey Zhang, Provided you have a no penetration boundary conditions for the airfoil, the general expression for the fluid force on an airfoil can be given as: F^i_R = \int_{\partial R} (P\hat{n} + \nu \nabla v^i) \cdot \hat{n} (where R is the airfoil domain, v^i is a unit vector in the i direction, google "Strict general mathematical definition of drag" for the complete derivation) Thanks. On Mon, Apr 27, 2015 at 4:15 AM, Zhang <zyzhang@...> wrote: Dear Vikram, Well. Just now I recognized that as for the form pressure force you are right. Since I want to get the lift and drag which are applied by the air flow onto the solid airfoil, the pressure form should be of positive. while in the original NS equation there is a negative one which means the force applied from solid to the gas, as opposite. In the same way, I think the contribution of viscosity onto the lift and drag should also be negative, since in the original NS form it is positive. Right? 原始邮件 发件人: "Vikram Garg" <vikram.v.garg@...> 发送时间: 20150427 16:32:10 (星期一) 收件人: Zhang <zyzhang@...> 抄送: libmeshusers <libmeshusers@...> 主题: Re: [Libmeshusers] calculation of aerodynamic forces and moments Hey Zhang, Yes, the get_normals() function will return the outward normal vectors. You had mentioned that force on the airfoil was given by: \vec{F}=\sum_{qpface_bound} {\tau*normal normal*p_{qp}} *Jxw_face_bound[qp] . Should that be a \tau*normal + normal*p_{qp}} ? Thanks. On Sun, Apr 26, 2015 at 11:07 PM, Zhang <zyzhang@...> wrote: I just made some correction, say, since pressure is defined in a space one order lower than velocity. In the last code I put all quadrature computation in velocity space by mistake, now I put the two integrals, $\int \tau * normal$ and $\int (p)*normal$, into two separate sections and finally add their partial integral together. Here are the comparison (for only pressure forces) 1) np=1 // the old one step 9, forces: (1.299297e01 3.460514e03 0.000000e+00 ) step 9, moments: (0.000000e+00 0.000000e+00 7.266930e05 ) 2) np=1 // the corrected one step 9, forces: (1.799026e01 3.255020e06 0.000000e+00 ) step 9, moments: (0.000000e+00 0.000000e+00 9.115135e07 ) 3) np=4 // the corrected one step 9, forces: (1.799026e01 2.098503e06 0.000000e+00 ) step 9, moments: (0.000000e+00 0.000000e+00 5.020619e07 ) when AOA=0.0 degree, the lift and zpitching moment are dramatically reduced to zero. Now the question is the negative xforce. In the face quadrature, does fe_pres_face>get_normals() return all outward normal vector? Zhenyu  One dashboard for servers and applications across PhysicalVirtualCloud Widest outofthebox monitoring support with 50+ applications Performance metrics, stats and reports that give you Actionable Insights Deep dive visibility with transaction tracing using APM Insight. http://ad.doubleclick.net/ddm/clk/290420510;117567292;y _______________________________________________ Libmeshusers mailing list Libmeshusers@... https://lists.sourceforge.net/lists/listinfo/libmeshusers  Vikram Garg Postdoctoral Associate Predictive Engineering and Computational Science (PECOS) The University of Texas at Austin http://web.mit.edu/vikramvg/www/ http://vikramvgarg.wordpress.com/ http://www.runforindia.org/runners/vikramg  One dashboard for servers and applications across PhysicalVirtualCloud Widest outofthebox monitoring support with 50+ applications Performance metrics, stats and reports that give you Actionable Insights Deep dive visibility with transaction tracing using APM Insight. http://ad.doubleclick.net/ddm/clk/290420510;117567292;y _______________________________________________ Libmeshusers mailing list Libmeshusers@... https://lists.sourceforge.net/lists/listinfo/libmeshusers  Vikram Garg Postdoctoral Associate Predictive Engineering and Computational Science (PECOS) The University of Texas at Austin http://web.mit.edu/vikramvg/www/ http://vikramvgarg.wordpress.com/ http://www.runforindia.org/runners/vikramg 
From: Vikram Garg <vikram.garg@gm...>  20150427 21:39:35

Hey Zhang, Provided you have a no penetration boundary conditions for the airfoil, the general expression for the fluid force on an airfoil can be given as: F^i_R = \int_{\partial R} (P\hat{n} + \nu \nabla v^i) \cdot \hat{n} (where R is the airfoil domain, v^i is a unit vector in the i direction, google "Strict general mathematical definition of drag" for the complete derivation) Thanks. On Mon, Apr 27, 2015 at 4:15 AM, Zhang <zyzhang@...> wrote: > Dear Vikram, > > Well. Just now I recognized that as for the form pressure force you are > right. > Since I want to get the lift and drag which are applied by the air flow > onto the solid airfoil, > the pressure form should be of positive. while in the original NS equation > there is a negative one which means the force applied from solid to the > gas, as opposite. > > In the same way, I think the contribution of viscosity onto the lift and > drag > should also be negative, since in the original NS form it is positive. > Right? > > 原始邮件 > 发件人: "Vikram Garg" <vikram.v.garg@...> > 发送时间: 20150427 16:32:10 (星期一) > 收件人: Zhang <zyzhang@...> > 抄送: libmeshusers <libmeshusers@...> > 主题: Re: [Libmeshusers] calculation of aerodynamic forces and moments > > > Hey Zhang, > Yes, the get_normals() function will return the outward > normal vectors. > > > You had mentioned that force on the airfoil was given by: > \vec{F}=\sum_{qpface_bound} {\tau*normal normal*p_{qp}} > *Jxw_face_bound[qp] . Should that be a \tau*normal + normal*p_{qp}} ? > > > Thanks. > > > On Sun, Apr 26, 2015 at 11:07 PM, Zhang <zyzhang@...> wrote: > I just made some correction, say, since pressure is defined in a space one > order > lower than velocity. In the last code I put all quadrature computation in > velocity space by mistake, > now I put the two integrals, $\int \tau * normal$ and $\int (p)*normal$, > into two separate > sections and finally add their partial integral together. > > Here are the comparison (for only pressure forces) > 1) np=1 // the old one > step 9, forces: (1.299297e01 3.460514e03 0.000000e+00 ) > step 9, moments: (0.000000e+00 0.000000e+00 7.266930e05 ) > 2) np=1 // the corrected one > step 9, forces: (1.799026e01 3.255020e06 0.000000e+00 ) > step 9, moments: (0.000000e+00 0.000000e+00 9.115135e07 ) > 3) np=4 // the corrected one > step 9, forces: (1.799026e01 2.098503e06 0.000000e+00 ) > step 9, moments: (0.000000e+00 0.000000e+00 5.020619e07 ) > > when AOA=0.0 degree, the lift and zpitching moment are dramatically > reduced to zero. > > Now the question is the negative xforce. In the face quadrature, does > fe_pres_face>get_normals() > return all outward normal vector? > > > Zhenyu > >  > One dashboard for servers and applications across PhysicalVirtualCloud > Widest outofthebox monitoring support with 50+ applications > Performance metrics, stats and reports that give you Actionable Insights > Deep dive visibility with transaction tracing using APM Insight. > http://ad.doubleclick.net/ddm/clk/290420510;117567292;y > _______________________________________________ > Libmeshusers mailing list > Libmeshusers@... > https://lists.sourceforge.net/lists/listinfo/libmeshusers > > > > > > >  > > Vikram Garg > Postdoctoral Associate > Predictive Engineering and Computational Science (PECOS) > The University of Texas at Austin > http://web.mit.edu/vikramvg/www/ > > > > http://vikramvgarg.wordpress.com/ > > http://www.runforindia.org/runners/vikramg > >  > One dashboard for servers and applications across PhysicalVirtualCloud > Widest outofthebox monitoring support with 50+ applications > Performance metrics, stats and reports that give you Actionable Insights > Deep dive visibility with transaction tracing using APM Insight. > http://ad.doubleclick.net/ddm/clk/290420510;117567292;y > _______________________________________________ > Libmeshusers mailing list > Libmeshusers@... > https://lists.sourceforge.net/lists/listinfo/libmeshusers >  Vikram Garg Postdoctoral Associate Predictive Engineering and Computational Science (PECOS) The University of Texas at Austin http://web.mit.edu/vikramvg/www/ http://vikramvgarg.wordpress.com/ http://www.runforindia.org/runners/vikramg 
From: Zhang <zyzhang@nu...>  20150427 09:31:37

Dear Vikram, Well. Just now I recognized that as for the form pressure force you are right. Since I want to get the lift and drag which are applied by the air flow onto the solid airfoil, the pressure form should be of positive. while in the original NS equation there is a negative one which means the force applied from solid to the gas, as opposite. In the same way, I think the contribution of viscosity onto the lift and drag should also be negative, since in the original NS form it is positive. Right? 原始邮件 发件人: "Vikram Garg" <vikram.v.garg@...> 发送时间: 20150427 16:32:10 (星期一) 收件人: Zhang <zyzhang@...> 抄送: libmeshusers <libmeshusers@...> 主题: Re: [Libmeshusers] calculation of aerodynamic forces and moments Hey Zhang, Yes, the get_normals() function will return the outward normal vectors. You had mentioned that force on the airfoil was given by: \vec{F}=\sum_{qpface_bound} {\tau*normal normal*p_{qp}} *Jxw_face_bound[qp] . Should that be a \tau*normal + normal*p_{qp}} ? Thanks. On Sun, Apr 26, 2015 at 11:07 PM, Zhang <zyzhang@...> wrote: I just made some correction, say, since pressure is defined in a space one order lower than velocity. In the last code I put all quadrature computation in velocity space by mistake, now I put the two integrals, $\int \tau * normal$ and $\int (p)*normal$, into two separate sections and finally add their partial integral together. Here are the comparison (for only pressure forces) 1) np=1 // the old one step 9, forces: (1.299297e01 3.460514e03 0.000000e+00 ) step 9, moments: (0.000000e+00 0.000000e+00 7.266930e05 ) 2) np=1 // the corrected one step 9, forces: (1.799026e01 3.255020e06 0.000000e+00 ) step 9, moments: (0.000000e+00 0.000000e+00 9.115135e07 ) 3) np=4 // the corrected one step 9, forces: (1.799026e01 2.098503e06 0.000000e+00 ) step 9, moments: (0.000000e+00 0.000000e+00 5.020619e07 ) when AOA=0.0 degree, the lift and zpitching moment are dramatically reduced to zero. Now the question is the negative xforce. In the face quadrature, does fe_pres_face>get_normals() return all outward normal vector? Zhenyu  One dashboard for servers and applications across PhysicalVirtualCloud Widest outofthebox monitoring support with 50+ applications Performance metrics, stats and reports that give you Actionable Insights Deep dive visibility with transaction tracing using APM Insight. http://ad.doubleclick.net/ddm/clk/290420510;117567292;y _______________________________________________ Libmeshusers mailing list Libmeshusers@... https://lists.sourceforge.net/lists/listinfo/libmeshusers  Vikram Garg Postdoctoral Associate Predictive Engineering and Computational Science (PECOS) The University of Texas at Austin http://web.mit.edu/vikramvg/www/ http://vikramvgarg.wordpress.com/ http://www.runforindia.org/runners/vikramg 
From: Vikram Garg <vikram.garg@gm...>  20150427 08:32:36

Hey Zhang, Yes, the get_normals() function will return the outward normal vectors. You had mentioned that force on the airfoil was given by: \vec{F}=\sum_{qpface_bound} {\tau*normal normal*p_{qp}} *Jxw_face_bound[qp] . Should that be a \tau*normal + normal*p_{qp}} ? Thanks. On Sun, Apr 26, 2015 at 11:07 PM, Zhang <zyzhang@...> wrote: > I just made some correction, say, since pressure is defined in a space one > order > lower than velocity. In the last code I put all quadrature computation in > velocity space by mistake, > now I put the two integrals, $\int \tau * normal$ and $\int (p)*normal$, > into two separate > sections and finally add their partial integral together. > > Here are the comparison (for only pressure forces) > 1) np=1 // the old one > step 9, forces: (1.299297e01 3.460514e03 0.000000e+00 ) > step 9, moments: (0.000000e+00 0.000000e+00 7.266930e05 ) > 2) np=1 // the corrected one > step 9, forces: (1.799026e01 3.255020e06 0.000000e+00 ) > step 9, moments: (0.000000e+00 0.000000e+00 9.115135e07 ) > 3) np=4 // the corrected one > step 9, forces: (1.799026e01 2.098503e06 0.000000e+00 ) > step 9, moments: (0.000000e+00 0.000000e+00 5.020619e07 ) > > when AOA=0.0 degree, the lift and zpitching moment are dramatically > reduced to zero. > > Now the question is the negative xforce. In the face quadrature, does > fe_pres_face>get_normals() > return all outward normal vector? > > Zhenyu > >  > One dashboard for servers and applications across PhysicalVirtualCloud > Widest outofthebox monitoring support with 50+ applications > Performance metrics, stats and reports that give you Actionable Insights > Deep dive visibility with transaction tracing using APM Insight. > http://ad.doubleclick.net/ddm/clk/290420510;117567292;y > _______________________________________________ > Libmeshusers mailing list > Libmeshusers@... > https://lists.sourceforge.net/lists/listinfo/libmeshusers >  Vikram Garg Postdoctoral Associate Predictive Engineering and Computational Science (PECOS) The University of Texas at Austin http://web.mit.edu/vikramvg/www/ http://vikramvgarg.wordpress.com/ http://www.runforindia.org/runners/vikramg 
From: Zhang <zyzhang@nu...>  20150427 04:22:20

I just made some correction, say, since pressure is defined in a space one order lower than velocity. In the last code I put all quadrature computation in velocity space by mistake, now I put the two integrals, $\int \tau * normal$ and $\int (p)*normal$, into two separate sections and finally add their partial integral together. Here are the comparison (for only pressure forces) 1) np=1 // the old one step 9, forces: (1.299297e01 3.460514e03 0.000000e+00 ) step 9, moments: (0.000000e+00 0.000000e+00 7.266930e05 ) 2) np=1 // the corrected one step 9, forces: (1.799026e01 3.255020e06 0.000000e+00 ) step 9, moments: (0.000000e+00 0.000000e+00 9.115135e07 ) 3) np=4 // the corrected one step 9, forces: (1.799026e01 2.098503e06 0.000000e+00 ) step 9, moments: (0.000000e+00 0.000000e+00 5.020619e07 ) when AOA=0.0 degree, the lift and zpitching moment are dramatically reduced to zero. Now the question is the negative xforce. In the face quadrature, does fe_pres_face>get_normals() return all outward normal vector? Zhenyu 
From: Zhang <zyzhang@nu...>  20150427 02:22:29

Dear Vikram, I just deal with product between tensors and vectors as a whole for these force/moment components together, not separating them before. If I properly understand your suggestion, do you think there's something wrong with the tensor computation? Anyway, I will tell you whether that's OK or not later for terms like \int{tau*normal }. Now, the results with only viscous shear stress are as follows NACA0012 airfoil, AOA=0., Re=68.8, 1) np=1 step 9, forces: (1.634845e+00 2.619208e07 0.000000e+00 ) step 9, moments: (0.000000e+00 0.000000e+00 8.167970e09 ) // Obviously the drag should NEVER be negative ! Lift and pitching moments are small enough yet. The full log is attached as result_np1.log. 2) np=4 step 9, forces: (1.634845e+00 1.288435e07 0.000000e+00 ) step 9, moments: (0.000000e+00 0.000000e+00 2.282703e08 ) still negative drag The full screen output are attached as from stdout.processor.0 to stdout.processor.3 here. In fact, my doubt comes from two points: 1) I may misused face quadrature for the boundary integral, esp. psi_face for p and q_face for velocity. could you please have a look at the code in my last email? 2) when running parallel with np 4, the big wrong result I got before (AOA>=4. ) always comes from a single node, rank 1 or 2, but this does not appear when AOA=0.0 with only viscous stress. Many thanks, Zhenyu 
From: Vikram Garg <vikram.garg@gm...>  20150427 00:15:23

Hello Zhang, What happens if you just try to compute the viscous drag (or vice versa) ? If one of the drags gives a reasonable answer and the other diverges, we can narrow things down. Thanks. On Sun, Apr 26, 2015 at 11:06 AM, Zhang <zyzhang@...> wrote: > Hi, > > I got some problems when writing some functions for normal, axial forces > and pitching moments of an airfoil in an incompressible laminar flow > (Re=68.8). > > Here I calculate the forces on the wing boundary elements as follows, > \vec{F}=\sum_{qpface_bound} {\tau*normal normal*p_{qp}} > *Jxw_face_bound[qp] > \vec{M}=\sum_{qpface_bound} R.cross({\tau*normal normal*p_{qp}}) > *Jxw_face_bound[qp] > where > p_{qp}=\sum_{0<=m<n_psi_face}{psi_face[m][qp]*soln(dof_p[m])}, > normal=(fe_vel_face>get_normals())[qp]; > tau=2*S*visc; > S=0.5*(gradu[m](n)+gradu[n](m)), m,n=0,1,2 > R=fe_vel_face>get_xyz()center_xyz; > > When I run the code, sometimes it gives quite large nonphysical force > components values like 4.698877e+256 , even at serial run. I guess I must > have made some improper usage of face quadratures. > Should I use only the velocity quadrature, or with the pressure one as > well? > > Here is the original code section for this work: > > template <unsigned int dimension> > void BoundaryMoment<dimension>::compute() > { > EquationSystems &es = *_equation_system; > > const MeshBase& mesh = es.get_mesh(); > > const unsigned int dim = mesh.mesh_dimension(); > > TransientNonlinearImplicitSystem& system = > es.get_system<TransientNonlinearImplicitSystem>(this>system_name); > > NumericVector<Number>& soln =*(system.current_local_solution); > const unsigned int u_var = system.variable_number ("u"); > const unsigned int p_var = system.variable_number ("p"); > > FEType fe_vel_type = system.variable_type(u_var); > FEType fe_pres_type = system.variable_type(p_var); > > const DofMap& dof_map = system.get_dof_map(); > > AutoPtr<FEBase> fe_vel_face (FEBase::build(dim, fe_vel_type)); > AutoPtr<FEBase> fe_pres_face (FEBase::build(dim, fe_pres_type)); > > QGauss qface(dim1, fe_vel_type.default_quadrature_order()); > QGauss qface_pres(dim1, fe_pres_type.default_quadrature_order()); > > fe_vel_face>attach_quadrature_rule (&qface); > fe_pres_face>attach_quadrature_rule (&qface_pres); > > const std::vector<std::vector<Real> >& phi_face = > fe_vel_face>get_phi(); > const std::vector<std::vector<RealGradient> >& dphi_face = > fe_vel_face>get_dphi(); > const std::vector<Real>& JxW_face = fe_vel_face>get_JxW(); > const std::vector<Point>& qface_normals = fe_vel_face>get_normals(); > const std::vector<Point>& qface_points = fe_vel_face>get_xyz(); > > const std::vector<std::vector<Real> >& psi_face = > fe_pres_face>get_phi(); > const std::vector<std::vector<RealGradient> >& dpsi_face = > fe_pres_face>get_dphi(); > const std::vector<Real>& JxW_face_pres = fe_pres_face>get_JxW(); > const std::vector<Point>& qface_normals_pres = > fe_pres_face>get_normals(); > const std::vector<Point>& qface_points_pres = fe_pres_face>get_xyz(); > > std::vector<dof_id_type> dof_indices; > std::vector<std::vector<dof_id_type> > dof_indices_u; > dof_indices_u.resize(this>_dim); > std::vector<dof_id_type> dof_indices_p; > > const Real visc = es.parameters.get<Real>("viscosity"); > > this>_values.zero(); > > const BoundaryInfo& b_info=mesh.get_boundary_info(); > > RealTensorValue _tau,_S2; > > Number p; > Number jxw; > > Point r; //vector from _origin to any qp point > > MeshBase::const_element_iterator el = > mesh.active_local_elements_begin(); > const MeshBase::const_element_iterator end_el = > mesh.active_local_elements_end(); > > int nelem=0; > for ( ; el != end_el; ++el) > { > const Elem* elem = *el; > > for(unsigned int d=0;d<this>_dim;d++) > dof_map.dof_indices (elem, dof_indices_u[d], u_var+d); > dof_map.dof_indices (elem, dof_indices_p, p_var); > > for (unsigned int side=0; side<elem>n_sides(); side++) > { > > bool boundary_found=false; > for(boundary_id_type& bid:this>_boundary_ids) > { > if(b_info.has_boundary_id(elem,side,bid)) > { > boundary_found=true; > break; > } > } > if(boundary_found) > { > AutoPtr<Elem> _side (elem>build_side(side)); > fe_vel_face>reinit(elem, side); > fe_pres_face>reinit(elem, side); > const std::vector<Point>& normals = > fe_vel_face>get_normals(); > > const std::vector<Point>& > qface_points=fe_vel_face>get_xyz(); // physical xyz > > for (unsigned int qp=0; qp<qface.n_points(); qp++) > { > Gradient gradu[3]; > r=qface_points[qp]_origin; > const Point& _normal=normals[qp]; > jxw=JxW_face[qp]; > > gradu[0].zero(); > gradu[1].zero(); > gradu[2].zero(); > for (unsigned int l=0; l<phi_face.size(); l++) > { > for(unsigned int d=0;d<this>_dim;d++) > gradu[d].add_scaled (dphi_face[l][qp],soln > (dof_indices_u[d][l])); > } > > _tau.zero(); > _S2.zero(); > > for(unsigned int m=0; m<this>_dim; m++) > for(unsigned int n=0; n<this>_dim; n++) > { > _S2(m,n)=gradu[m](n); > } > _S2+=_S2.transpose(); > _tau=_S2*visc; > > p=0.; > // for (unsigned int l=0; l<n_p_dofs; l++) > for (unsigned int l=0; l<psi_face.size(); l++)//??? > { > p += psi_face[l][qp]*soln(dof_indices_p[l]); > } > > this>_values+=r.cross(_tau*_normal_normal*p)*jxw; > }// qp loop in qface > nelem++; > } // if (elem>on_boundary()) > }// for (unsigned int side=0; side<elem>n_sides(); side++) > } // el loop > > #if 1 > system.comm().sum<Number>(this>_values(0)); > system.comm().sum<Number>(this>_values(1)); > system.comm().sum<Number>(this>_values(2)); > #endif > > cout<<"nelem=: "<<nelem<<"; _values: " > <<this>_values(0)<<"," > <<this>_values(1)<<"," > <<this>_values(2) > <<endl; > > } // void BoundaryMoment::compute() > > Thank you for any possible help, > > Zhenyu > >  > One dashboard for servers and applications across PhysicalVirtualCloud > Widest outofthebox monitoring support with 50+ applications > Performance metrics, stats and reports that give you Actionable Insights > Deep dive visibility with transaction tracing using APM Insight. > http://ad.doubleclick.net/ddm/clk/290420510;117567292;y > _______________________________________________ > Libmeshusers mailing list > Libmeshusers@... > https://lists.sourceforge.net/lists/listinfo/libmeshusers >  Vikram Garg Postdoctoral Associate Predictive Engineering and Computational Science (PECOS) The University of Texas at Austin http://web.mit.edu/vikramvg/www/ http://vikramvgarg.wordpress.com/ http://www.runforindia.org/runners/vikramg 
From: Zhang <zyzhang@nu...>  20150426 16:22:15

Hi, I got some problems when writing some functions for normal, axial forces and pitching moments of an airfoil in an incompressible laminar flow (Re=68.8). Here I calculate the forces on the wing boundary elements as follows, \vec{F}=\sum_{qpface_bound} {\tau*normal normal*p_{qp}} *Jxw_face_bound[qp] \vec{M}=\sum_{qpface_bound} R.cross({\tau*normal normal*p_{qp}}) *Jxw_face_bound[qp] where p_{qp}=\sum_{0<=m<n_psi_face}{psi_face[m][qp]*soln(dof_p[m])}, normal=(fe_vel_face>get_normals())[qp]; tau=2*S*visc; S=0.5*(gradu[m](n)+gradu[n](m)), m,n=0,1,2 R=fe_vel_face>get_xyz()center_xyz; When I run the code, sometimes it gives quite large nonphysical force components values like 4.698877e+256 , even at serial run. I guess I must have made some improper usage of face quadratures. Should I use only the velocity quadrature, or with the pressure one as well? Here is the original code section for this work: template <unsigned int dimension> void BoundaryMoment<dimension>::compute() { EquationSystems &es = *_equation_system; const MeshBase& mesh = es.get_mesh(); const unsigned int dim = mesh.mesh_dimension(); TransientNonlinearImplicitSystem& system = es.get_system<TransientNonlinearImplicitSystem>(this>system_name); NumericVector<Number>& soln =*(system.current_local_solution); const unsigned int u_var = system.variable_number ("u"); const unsigned int p_var = system.variable_number ("p"); FEType fe_vel_type = system.variable_type(u_var); FEType fe_pres_type = system.variable_type(p_var); const DofMap& dof_map = system.get_dof_map(); AutoPtr<FEBase> fe_vel_face (FEBase::build(dim, fe_vel_type)); AutoPtr<FEBase> fe_pres_face (FEBase::build(dim, fe_pres_type)); QGauss qface(dim1, fe_vel_type.default_quadrature_order()); QGauss qface_pres(dim1, fe_pres_type.default_quadrature_order()); fe_vel_face>attach_quadrature_rule (&qface); fe_pres_face>attach_quadrature_rule (&qface_pres); const std::vector<std::vector<Real> >& phi_face = fe_vel_face>get_phi(); const std::vector<std::vector<RealGradient> >& dphi_face = fe_vel_face>get_dphi(); const std::vector<Real>& JxW_face = fe_vel_face>get_JxW(); const std::vector<Point>& qface_normals = fe_vel_face>get_normals(); const std::vector<Point>& qface_points = fe_vel_face>get_xyz(); const std::vector<std::vector<Real> >& psi_face = fe_pres_face>get_phi(); const std::vector<std::vector<RealGradient> >& dpsi_face = fe_pres_face>get_dphi(); const std::vector<Real>& JxW_face_pres = fe_pres_face>get_JxW(); const std::vector<Point>& qface_normals_pres = fe_pres_face>get_normals(); const std::vector<Point>& qface_points_pres = fe_pres_face>get_xyz(); std::vector<dof_id_type> dof_indices; std::vector<std::vector<dof_id_type> > dof_indices_u; dof_indices_u.resize(this>_dim); std::vector<dof_id_type> dof_indices_p; const Real visc = es.parameters.get<Real>("viscosity"); this>_values.zero(); const BoundaryInfo& b_info=mesh.get_boundary_info(); RealTensorValue _tau,_S2; Number p; Number jxw; Point r; //vector from _origin to any qp point MeshBase::const_element_iterator el = mesh.active_local_elements_begin(); const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end(); int nelem=0; for ( ; el != end_el; ++el) { const Elem* elem = *el; for(unsigned int d=0;d<this>_dim;d++) dof_map.dof_indices (elem, dof_indices_u[d], u_var+d); dof_map.dof_indices (elem, dof_indices_p, p_var); for (unsigned int side=0; side<elem>n_sides(); side++) { bool boundary_found=false; for(boundary_id_type& bid:this>_boundary_ids) { if(b_info.has_boundary_id(elem,side,bid)) { boundary_found=true; break; } } if(boundary_found) { AutoPtr<Elem> _side (elem>build_side(side)); fe_vel_face>reinit(elem, side); fe_pres_face>reinit(elem, side); const std::vector<Point>& normals = fe_vel_face>get_normals(); const std::vector<Point>& qface_points=fe_vel_face>get_xyz(); // physical xyz for (unsigned int qp=0; qp<qface.n_points(); qp++) { Gradient gradu[3]; r=qface_points[qp]_origin; const Point& _normal=normals[qp]; jxw=JxW_face[qp]; gradu[0].zero(); gradu[1].zero(); gradu[2].zero(); for (unsigned int l=0; l<phi_face.size(); l++) { for(unsigned int d=0;d<this>_dim;d++) gradu[d].add_scaled (dphi_face[l][qp],soln (dof_indices_u[d][l])); } _tau.zero(); _S2.zero(); for(unsigned int m=0; m<this>_dim; m++) for(unsigned int n=0; n<this>_dim; n++) { _S2(m,n)=gradu[m](n); } _S2+=_S2.transpose(); _tau=_S2*visc; p=0.; // for (unsigned int l=0; l<n_p_dofs; l++) for (unsigned int l=0; l<psi_face.size(); l++)//??? { p += psi_face[l][qp]*soln(dof_indices_p[l]); } this>_values+=r.cross(_tau*_normal_normal*p)*jxw; }// qp loop in qface nelem++; } // if (elem>on_boundary()) }// for (unsigned int side=0; side<elem>n_sides(); side++) } // el loop #if 1 system.comm().sum<Number>(this>_values(0)); system.comm().sum<Number>(this>_values(1)); system.comm().sum<Number>(this>_values(2)); #endif cout<<"nelem=: "<<nelem<<"; _values: " <<this>_values(0)<<"," <<this>_values(1)<<"," <<this>_values(2) <<endl; } // void BoundaryMoment::compute() Thank you for any possible help, Zhenyu 
From: Xikai Jiang <jxk28@ms...>  20150424 20:01:20

Thanks for your prompt response. I just realize the last post about "infinite element" on the mailing list is a decade ago.. It looks like Steffen were at Standford at that time. I'll go through the webpages you point out, and see what I can do. Xikai From: jwpeterson@... Date: Fri, 24 Apr 2015 13:44:08 0600 Subject: Re: [Libmeshusers] infinite element Laplace To: ptbauman@... CC: jxk28@...; libmeshusers@... On Fri, Apr 24, 2015 at 1:42 PM, Paul T. Bauman <ptbauman@...> wrote: On Fri, Apr 24, 2015 at 3:34 PM, Xikai Jiang <jxk28@...> wrote: > > 1. Could I use the "mapped wave envelop element" for 3D Laplace equation > in infinite domain? i.e. to accurately describe 1/r decay? > > 2. If yes, what would you suggest in order to use "mapped wave envelop > element" for 3D Laplace? > > 3. If no, are there any types of infinite element implemented in libMesh > for 3D Laplace? > Others will hopefully know more than me, but there are 3D infinite elements, both geometric http://libmesh.github.io/doxygen/classlibMesh_1_1InfCell.html and finite elements http://libmesh.github.io/doxygen/classlibMesh_1_1InfFE.html I have no idea how they work. The person that wrote those has not been active in libMesh development for a long while. The original infinite element developers were Daniel Dreyer and Steffen Petersen, but I currently don't know how to get in touch with either of them.  John 
From: John Peterson <jwpeterson@gm...>  20150424 19:44:34

On Fri, Apr 24, 2015 at 1:42 PM, Paul T. Bauman <ptbauman@...> wrote: > On Fri, Apr 24, 2015 at 3:34 PM, Xikai Jiang <jxk28@...> wrote: > > > > > 1. Could I use the "mapped wave envelop element" for 3D Laplace equation > > in infinite domain? i.e. to accurately describe 1/r decay? > > > > 2. If yes, what would you suggest in order to use "mapped wave envelop > > element" for 3D Laplace? > > > > 3. If no, are there any types of infinite element implemented in libMesh > > for 3D Laplace? > > > > Others will hopefully know more than me, but there are 3D infinite > elements, both geometric > http://libmesh.github.io/doxygen/classlibMesh_1_1InfCell.html and finite > elements http://libmesh.github.io/doxygen/classlibMesh_1_1InfFE.html > > I have no idea how they work. The person that wrote those has not been > active in libMesh development for a long while. > The original infinite element developers were Daniel Dreyer and Steffen Petersen, but I currently don't know how to get in touch with either of them.  John 
From: Paul T. Bauman <ptbauman@gm...>  20150424 19:42:32

On Fri, Apr 24, 2015 at 3:34 PM, Xikai Jiang <jxk28@...> wrote: > > 1. Could I use the "mapped wave envelop element" for 3D Laplace equation > in infinite domain? i.e. to accurately describe 1/r decay? > > 2. If yes, what would you suggest in order to use "mapped wave envelop > element" for 3D Laplace? > > 3. If no, are there any types of infinite element implemented in libMesh > for 3D Laplace? > Others will hopefully know more than me, but there are 3D infinite elements, both geometric http://libmesh.github.io/doxygen/classlibMesh_1_1InfCell.html and finite elements http://libmesh.github.io/doxygen/classlibMesh_1_1InfFE.html I have no idea how they work. The person that wrote those has not been active in libMesh development for a long while. I'm hopeful some of the other developers have more knowledge. I can at least point you to an example where they're used (wave equation): http://libmesh.github.io/examples/miscellaneous_ex1.html 
From: Xikai Jiang <jxk28@ms...>  20150424 19:34:26

Dear libMesh Users and Developers, I want to simulate a 3D electrostatic problem in infinite domain. The electrical potential follows 3D Laplace equation, and it decays to zero (in the form of 1/r) or reaches a finite value at infinity. I'm looking into infinite element in libMesh. By searching the mailing list archive, I got to know the infinite element implemented in libMesh is "mapped wave envelop element" for acoustic radiation and scattering, my questions are, 1. Could I use the "mapped wave envelop element" for 3D Laplace equation in infinite domain? i.e. to accurately describe 1/r decay? 2. If yes, what would you suggest in order to use "mapped wave envelop element" for 3D Laplace? 3. If no, are there any types of infinite element implemented in libMesh for 3D Laplace? Thanks for your time. Regards, Xikai 
From: Paul T. Bauman <ptbauman@gm...>  20150424 12:56:16

replyall fail On Fri, Apr 24, 2015 at 8:54 AM, David Knezevic <david.knezevic@...> wrote: > > Sounds good to me. I can make a PR for this if you'd like? > That would be great, thanks! best, Paul 
From: David Knezevic <david.knezevic@ak...>  20150424 12:54:48

> > Adding add _elem_dims.clear(); at the start of MeshBase::cache_elem_dims() >> seems to fix the problem for me. >> > > Yeah, that's the fix. We should probably add some documentation to > delete_elem() (and add_elem() for that matter) that the user should call > cache_elem_dims() when they're done (as you did in your code). > Sounds good to me. I can make a PR for this if you'd like? Or feel free to make the change yourself if you prefer. David 
From: Paul T. Bauman <ptbauman@gm...>  20150424 12:37:54

On Fri, Apr 24, 2015 at 8:28 AM, David Knezevic <david.knezevic@...> wrote: > > Not sure if you saw my followup email? > Yeah, sorry, I just wanted to make sure what the code path you were following was before speculating. > Adding add _elem_dims.clear(); at the start of MeshBase::cache_elem_dims() > seems to fix the problem for me. > Yeah, that's the fix. We should probably add some documentation to delete_elem() (and add_elem() for that matter) that the user should call cache_elem_dims() when they're done (as you did in your code). 
From: Paul T. Bauman <ptbauman@gm...>  20150424 12:29:21

On Fri, Apr 24, 2015 at 12:00 AM, David Knezevic <david.knezevic@... > wrote: > I've attached code that does this, with an associated test mesh. Could send to me (and Roy for good measure since he's participating in the thread) the repro case? Looks like the list stripped it. Thanks! 
From: David Knezevic <david.knezevic@ak...>  20150424 12:28:47

On Fri, Apr 24, 2015 at 6:39 AM, Paul T. Bauman <ptbauman@...> wrote: > > > > On Apr 24, 2015, at 12:00 AM, David Knezevic <david.knezevic@...> > wrote: > > > > I'd like to be able to impose boundary conditions on an edge. We added > > support for this in libMesh, but the only way I could come up with for > > creating edge BCs in a mesh generator is to: > > Two ideas, neither of which I've tested: > > 1. You ought to be able to select curves in Cubit (IIRC you were using > this one?) when building side sets. You probably already tried this, but > wanted to be doubly sure. Does it dump them from the sideset because it's > not a face? > As far as we could tell, this doesn't work in CUBIT. I asked someone from CUBIT about it and he confirmed that CUBIT doesn't support sidesets for edges. > 2. I've looked at the code in libmesh, but again haven't tested, but if > you only need nodal dofs constrained, you could create a node set and treat > it like a boundary id. The libmesh code looked like it effectively treated > node sets like side sets. Maybe you could use this too to find nonnodal > dofs on the edge and add them to the node set once it's been created? > This could work as a fallback, but it seems like building libMesh edge BCs is nicer since that automatically handles any edge dofs. The purpose of those 1D elements in the mesh is just to allow us to construct edge BCs in the correct place. I'll also have a quick look later this morning at the bug you found. Not sure if you saw my followup email? Adding add _elem_dims.clear(); at the start of MeshBase::cache_elem_dims() seems to fix the problem for me. David 