You can subscribe to this list here.
2003 
_{Jan}
(4) 
_{Feb}
(1) 
_{Mar}
(9) 
_{Apr}
(2) 
_{May}
(7) 
_{Jun}
(1) 
_{Jul}
(1) 
_{Aug}
(4) 
_{Sep}
(12) 
_{Oct}
(8) 
_{Nov}
(3) 
_{Dec}
(4) 

2004 
_{Jan}
(1) 
_{Feb}
(21) 
_{Mar}
(31) 
_{Apr}
(10) 
_{May}
(12) 
_{Jun}
(15) 
_{Jul}
(4) 
_{Aug}
(6) 
_{Sep}
(5) 
_{Oct}
(11) 
_{Nov}
(43) 
_{Dec}
(13) 
2005 
_{Jan}
(25) 
_{Feb}
(12) 
_{Mar}
(49) 
_{Apr}
(19) 
_{May}
(104) 
_{Jun}
(60) 
_{Jul}
(10) 
_{Aug}
(42) 
_{Sep}
(15) 
_{Oct}
(12) 
_{Nov}
(6) 
_{Dec}
(4) 
2006 
_{Jan}
(1) 
_{Feb}
(6) 
_{Mar}
(31) 
_{Apr}
(17) 
_{May}
(5) 
_{Jun}
(95) 
_{Jul}
(38) 
_{Aug}
(44) 
_{Sep}
(6) 
_{Oct}
(8) 
_{Nov}
(21) 
_{Dec}

2007 
_{Jan}
(5) 
_{Feb}
(46) 
_{Mar}
(9) 
_{Apr}
(23) 
_{May}
(17) 
_{Jun}
(51) 
_{Jul}
(41) 
_{Aug}
(4) 
_{Sep}
(28) 
_{Oct}
(71) 
_{Nov}
(193) 
_{Dec}
(20) 
2008 
_{Jan}
(46) 
_{Feb}
(46) 
_{Mar}
(18) 
_{Apr}
(38) 
_{May}
(14) 
_{Jun}
(107) 
_{Jul}
(50) 
_{Aug}
(115) 
_{Sep}
(84) 
_{Oct}
(96) 
_{Nov}
(105) 
_{Dec}
(34) 
2009 
_{Jan}
(89) 
_{Feb}
(93) 
_{Mar}
(119) 
_{Apr}
(73) 
_{May}
(39) 
_{Jun}
(51) 
_{Jul}
(27) 
_{Aug}
(8) 
_{Sep}
(91) 
_{Oct}
(90) 
_{Nov}
(77) 
_{Dec}
(67) 
2010 
_{Jan}
(25) 
_{Feb}
(36) 
_{Mar}
(98) 
_{Apr}
(45) 
_{May}
(25) 
_{Jun}
(60) 
_{Jul}
(17) 
_{Aug}
(36) 
_{Sep}
(48) 
_{Oct}
(45) 
_{Nov}
(65) 
_{Dec}
(39) 
2011 
_{Jan}
(26) 
_{Feb}
(48) 
_{Mar}
(151) 
_{Apr}
(108) 
_{May}
(61) 
_{Jun}
(108) 
_{Jul}
(27) 
_{Aug}
(50) 
_{Sep}
(43) 
_{Oct}
(43) 
_{Nov}
(27) 
_{Dec}
(37) 
2012 
_{Jan}
(56) 
_{Feb}
(120) 
_{Mar}
(72) 
_{Apr}
(57) 
_{May}
(82) 
_{Jun}
(66) 
_{Jul}
(51) 
_{Aug}
(75) 
_{Sep}
(166) 
_{Oct}
(232) 
_{Nov}
(284) 
_{Dec}
(105) 
2013 
_{Jan}
(168) 
_{Feb}
(151) 
_{Mar}
(30) 
_{Apr}
(145) 
_{May}
(26) 
_{Jun}
(53) 
_{Jul}
(76) 
_{Aug}
(33) 
_{Sep}
(23) 
_{Oct}
(72) 
_{Nov}
(125) 
_{Dec}
(38) 
2014 
_{Jan}
(47) 
_{Feb}
(62) 
_{Mar}
(27) 
_{Apr}
(8) 
_{May}
(12) 
_{Jun}
(2) 
_{Jul}
(22) 
_{Aug}
(22) 
_{Sep}

_{Oct}
(17) 
_{Nov}
(20) 
_{Dec}
(12) 
2015 
_{Jan}
(25) 
_{Feb}
(2) 
_{Mar}
(16) 
_{Apr}
(13) 
_{May}
(21) 
_{Jun}
(5) 
_{Jul}
(1) 
_{Aug}
(8) 
_{Sep}
(9) 
_{Oct}
(30) 
_{Nov}
(8) 
_{Dec}

2016 
_{Jan}
(16) 
_{Feb}
(31) 
_{Mar}
(43) 
_{Apr}
(18) 
_{May}
(21) 
_{Jun}
(11) 
_{Jul}
(17) 
_{Aug}
(26) 
_{Sep}

_{Oct}

_{Nov}

_{Dec}

S  M  T  W  T  F  S 

1

2

3
(1) 
4

5
(1) 
6

7

8
(4) 
9
(6) 
10
(3) 
11
(1) 
12
(1) 
13

14

15

16

17

18

19

20

21

22

23

24

25

26

27
(4) 
28

29







From: Steffen Petersen <steffen.petersen@tu...>  20040227 20:17:12

With some limitations I got libmesh working on our HP machines :). Since log seems to be already defined in the hpuxstl I would like to change libMesh::log to e.g. libMesh::logging. Any objections? Steffen 
From: William L. (Bill) Barth <bbarth@cf...>  20040227 16:04:25

>>>>> On Fri, 27 Feb 2004 08:59:33 0600, "KIRK, BENJAMIN (JSCEG) (NASA)" <benjamin.kirk1@...> said: Ben> So foo1 should probably be called normal(), but what about foo2? I Ben> suggested size() some time back, but as John pointed out that name generally Ben> has a different meaning thanks to the STL. What has a similar meaning to Ben> length/area/volume but is dimensionindependent?? I'd just call it volume() or vol() since it's atop the language hierarchy, and then put a comment in the documentation along the lines of: Returns the volume/area/length of side s or the element itself with no argument. Bill.  Bill Barth  Home: (512) 7973045 bbarth@...  Work: (512) 4714069 Office: WRW 111  Fax: (512) 2323357 
From: John Peterson <peterson@cf...>  20040227 15:56:12

KIRK, BENJAMIN (JSCEG) (NASA) writes: > I need some help thinking of a good method name... > > Here is what I want: > > Point normal = elem>foo1(s); > // returns a vector normal to side s. > // the magnitude of this vector is the area > // of side s. > > Real lav = elem>foo2(s); > // returns the length or area of side s (2D, 3D respectively) > // s defaults to libMesh::invalid_uint, in which case > // elem>foo2() returns the length/area/volume of the element. > > So foo1 should probably be called normal(), but what about foo2? I > suggested size() some time back, but as John pointed out that name generally > has a different meaning thanks to the STL. What has a similar meaning to > length/area/volume but is dimensionindependent?? foo2() is dimension independent ;) but I would suggest elem>size_of_side(s) It's longer but tells exactly what it's doing... John 
From: KIRK, BENJAMIN (JSCEG) (NASA) <benjamin.kirk1@na...>  20040227 15:10:01

I need some help thinking of a good method name... Here is what I want: Point normal = elem>foo1(s); // returns a vector normal to side s. // the magnitude of this vector is the area // of side s. Real lav = elem>foo2(s); // returns the length or area of side s (2D, 3D respectively) // s defaults to libMesh::invalid_uint, in which case // elem>foo2() returns the length/area/volume of the element. So foo1 should probably be called normal(), but what about foo2? I suggested size() some time back, but as John pointed out that name generally has a different meaning thanks to the STL. What has a similar meaning to length/area/volume but is dimensionindependent?? Ben 
From: KIRK, BENJAMIN (JSCEG) (NASA) <benjamin.kirk1@na...>  20040212 17:42:28

Does anyone have any presentations or publications to add to the website? They don't have to reference libMesh explicitly... I would just like to get some more stuff up there. If the issue is editing the PHP so that the page works right then don't worry... Either John or I can do that if you don't want to. Ben 
From: Steffen Petersen <steffen.petersen@tu...>  20040211 10:57:22

I can give it a shot on a HP machine, but it is quite some time ago that Daniel got the library running on our HPs. Steffen > I'm testing portability now, but I don't have access to every > machine we claim to support (particularly HP). Will anyone > with access to any "problem" machine (where libMesh *has* run > in the past, but you haven't tried recently) try to build the > library? Don't worry about GNULinux machines... > > Thanks. > > Ben 
From: William L. (Bill) Barth <bbarth@cf...>  20040210 20:45:53

Ben, can you give a status update on the planned full decomposition/parallelization of the mesh? Bill.  Bill Barth  Home: (512) 7973045 bbarth@...  Work: (512) 4714069 Office: WRW 111  Fax: (512) 2323357 
From: KIRK, BENJAMIN (JSCEG) (NASA) <benjamin.kirk1@na...>  20040210 20:07:46

I'm testing portability now, but I don't have access to every machine we claim to support (particularly HP). Will anyone with access to any "problem" machine (where libMesh *has* run in the past, but you haven't tried recently) try to build the library? Don't worry about GNULinux machines... Thanks. Ben 
From: Steffen Petersen <steffen.petersen@tu...>  20040210 16:55:39

I have been away for a few days so I did not get any mail until today. I just commited szabab shapes with nicer representations (similar to what I did in 1D). Hendrik had already converted them last week but we did not commit the changes since they were not tested. Some pow<2>s are still in the code which I think should be ok. So the szabab elements should now be ready for the new release :) Steffen > Daniel Dreyer writes: > > I've seen you already changed the sources  Thanks a lot! > > > > Please help me then to clarify something: i think i picked > it > up in the deal.ii library, that Horner is best w.r.t. > > avoiding roundoff, simply because 0.1^5 = 1.e5, and > > when you may substract this from something else of order > 1.e1, > you may get orderofmagnitude cancellation. That is > true. I see that Steffen (or somebody) already implemented > the Horner scheme for the 1D routines. If anyone wants to do > this for the 2D shapes I highly encourage (and pity!) them. > By putting all the pow<> in there, I was not trying to imply > that was the "correct" way to do things :) > > John > > >  > The SF.Net email is sponsored by EclipseCon 2004 > Premiere Conference on Open Tools Development and Integration > See the breadth of Eclipse activity. February 35 in Anaheim, > CA. http://www.eclipsecon.org/osdn > _______________________________________________ > Libmeshdevel mailing list > Libmeshdevel@... > https://lists.sourceforge.net/lists/listinfo/libmeshdevel > 
From: John Peterson <peterson@cf...>  20040209 20:35:47

Daniel Dreyer writes: > I've seen you already changed the sources  Thanks a lot! > > Please help me then to clarify something: i think i picked it > up in the deal.ii library, that Horner is best w.r.t. > avoiding roundoff, simply because 0.1^5 = 1.e5, and > when you may substract this from something else of order 1.e1, > you may get orderofmagnitude cancellation. That is true. I see that Steffen (or somebody) already implemented the Horner scheme for the 1D routines. If anyone wants to do this for the 2D shapes I highly encourage (and pity!) them. By putting all the pow<> in there, I was not trying to imply that was the "correct" way to do things :) John 
From: KIRK, BENJAMIN (JSCEG) (NASA) <benjamin.kirk1@na...>  20040209 18:13:54

As I remember it, Hughes, Tezduyar, and others use an elementbyelement matrixfree GMRES approach. In this approach the global matrix is never formed, and each time a matrixvector product is needed the element matrix assembly routine (or a functional equivalent) is called. This is more memoryefficient but more computationallyexpensive that assembling the global matrix and preconditioning it. In libMesh the global matrix is assembled, and it may be preconditioned as you see fit. This is generally better than elementbyelement preconditioning. However, if you want to perform elementbyelement preconditioning you can certainly precondition the element matrix (after you assemble it in your code) before you insert it into the global matrix. As for matrixfree, Daniel has expressed an interest in this in the past. It would require some extensions to the linear solver interface, but that is doable. Does this help? Ben Original Message From: libmeshdeveladmin@... [mailto:libmeshdeveladmin@...] On Behalf Of seid mehdi bostandoust nik Sent: Monday, February 09, 2004 12:02 PM To: libmeshusers@...; libmeshdevel@... Subject: [Libmeshdevel] element by element preconditioning? Hi Based on hughes,Tezduyar and others it is possible to solve NavierStokes equations with iterative solvers by the help of element by element preconditioning. I am using petsc and currently in Petsc there is not element by element preconditioning(it is said that this feature may be added to PETSC3.0). I want to know that does the libmesh solver support element by element preconditioning for iterative solvers or not? thanks bye __________________________________ Do you Yahoo!? Yahoo! Finance: Get your refund fast by filing online. http://taxes.yahoo.com/filing.html  The SF.Net email is sponsored by EclipseCon 2004 Premiere Conference on Open Tools Development and Integration See the breadth of Eclipse activity. February 35 in Anaheim, CA. http://www.eclipsecon.org/osdn _______________________________________________ Libmeshdevel mailing list Libmeshdevel@... https://lists.sourceforge.net/lists/listinfo/libmeshdevel 
From: seid mehdi bostandoust nik <mbostandoust@ya...>  20040209 18:01:46

Hi Based on hughes,Tezduyar and others it is possible to solve NavierStokes equations with iterative solvers by the help of element by element preconditioning. I am using petsc and currently in Petsc there is not element by element preconditioning(it is said that this feature may be added to PETSC3.0). I want to know that does the libmesh solver support element by element preconditioning for iterative solvers or not? thanks bye __________________________________ Do you Yahoo!? Yahoo! Finance: Get your refund fast by filing online. http://taxes.yahoo.com/filing.html 
From: John Peterson <peterson@cf...>  20040209 15:19:07

Benjamin S. Kirk writes: > Yup... that seems about right... > > The prototype for pow looks like: > > double pow(double, double); > > So there is no distinction between x^10 and x^pi, for instance... (Or if > there is a distinction it is implemented in C with iftests, which is > slow too. As John noted, we can do better assuming double^int with a > recursive template that is unrolled at compiletime. We can do even > better with a recursive template that handles repeated squaring, but for > small powers (< ~10) I don't expect to see much difference... > > The only change I might make to the code below is a provision to handle > negative powers... Maybe 2 templaes, one <int > that handles +/ > powers, and one <unsigned int> optimized for positive powers? > For negative powers just return the inverse of the positive power... One additional note if you are trying this at home. Add the inline keyword to both recursive template functions. It will not be nearly as fast if your compiler chooses not to inline those functions. John 
From: Benjamin S. Kirk <benkirk@cf...>  20040209 13:44:38

Yup... that seems about right... The prototype for pow looks like: double pow(double, double); So there is no distinction between x^10 and x^pi, for instance... (Or if there is a distinction it is implemented in C with iftests, which is slow too. As John noted, we can do better assuming double^int with a recursive template that is unrolled at compiletime. We can do even better with a recursive template that handles repeated squaring, but for small powers (< ~10) I don't expect to see much difference... The only change I might make to the code below is a provision to handle negative powers... Maybe 2 templaes, one <int > that handles +/ powers, and one <unsigned int> optimized for positive powers? For negative powers just return the inverse of the positive power... Ben John Peterson wrote: >Benjamin S. Kirk writes: > > Help me out here, John... Do you remeber that C++ User's Journal on > > scientific programming? I think in near the end of that they describe a > > templated pow function that hadles x^int cases intelligently through > > repeated squaring... Does that ring a bell? It might be worth having > > in utility.h > >OK, I've done some testing. Apparently pow() is ridiculously slow. >I tested two different methods and pow in the attached code. Under >O3 optimization, I found the results on a P4 2GHz (in seconds): > >Raising 50 million double values to the third power. >Recursive Template: 0.05 >Template Loop: 1.93 >Standard pow(): 13.01 > >Raising to the fifth power, results were: >Recursive Template: 0.05 >Template Loop: 2.23 >Standard pow(): 14.17 > >Finally, raising to the 10th power (base numbers were small) >Recursive Template: 0.06 >Template Loop: 2.64 >Standard pow(): 13.56 > >The "Recursive Template" method is the method I posted to the mailing >list earlier today. It seems to be the fastest, and as expected, its >performance remains about constant no matter how high the exponent becomes. >Here is the code (don't forget lm for pow()) I invite you to test it >out yourselves... > >#include <iostream> >#include <vector> >#include <time.h> >#include <iomanip> > > >// Anonymous Namespace to hold useful constants >namespace { >const unsigned int n_bases = static_cast<unsigned int>(1e7); > >// Number of tests to perform >const unsigned int n_tests = 5; > >// Raise to an integer power >const unsigned int power = 10; >}; > > > >// Power raising class >template <int N> >struct Pow >{ > template <int Val> > struct In > { > enum { value = Val * Pow<N1>::template In<Val>::value }; > }; >}; > >template <> >struct Pow<0> >{ > template <int Val> > struct In > { > enum { value = 1 }; > }; >}; > > > >// Rescursive template function >template <int N> >double pow(double x) { return x * pow<N1>(x); } >template <> >double pow<0>(double x) { return 1.; } > > >// Templated Loop function >template <int N> >double pow2(double x) >{ > const unsigned int l = N / 2; > const unsigned int m = N % 2; > const double x2 = x*x; > double y = 1.0; > > for (unsigned int i=0; i<l; ++i) > { > y = y * x2; > } > > return m==0 ? y : y*x; >} > >// I can't figure out how to get the class to use >// double as its base value. Therefore, this method >// is not feasible now. >// void test_class(std::vector<double>& base_values) >// { >// for (unsigned int n=0; n<n_tests; ++n) >// { >// for (unsigned int i=0; i<n_bases; ++i) >// //std::cout << >// Pow<base_values[i]>::In<power> >// // << std::endl >// ; >// } >// } > >void test_recursive_template(std::vector<double>& base_values) >{ > for (unsigned int n=0; n<n_tests; ++n) > { > for (unsigned int i=0; i<n_bases; ++i) > //std::cout << > pow<power>(base_values[i]) > // << std::endl > ; > } > >} > >void test_template_loop(std::vector<double>& base_values) >{ > for (unsigned int n=0; n<n_tests; ++n) > { > for (unsigned int i=0; i<n_bases; ++i) > //std::cout << > pow2<power>(base_values[i]) > // << std::endl > ; > } >} > > > >void test_standard_pow(std::vector<double>& base_values) >{ > for (unsigned int n=0; n<n_tests; ++n) > { > for (unsigned int i=0; i<n_bases; ++i) > //std::cout << > pow(base_values[i], power) > // << std::endl > ; > } >} > > > > > >int main() >{ > srand(0); > > std::vector<double> base_values(n_bases); > > // Fill the vector with random values between zero and one. > for (unsigned int i=0; i<n_bases; ++i) > base_values[i] = static_cast<double>(rand()) / > static_cast<double>(RAND_MAX); > > // Convenient typedef for function pointers > typedef void (*TestFunction) (std::vector<double>& v); > > // A vector of function pointers for the functions we will test. > std::vector<TestFunction> tests(3); > tests[0] = test_recursive_template; > tests[1] = test_template_loop; > tests[2] = test_standard_pow; > > // A vector of convenient names > std::vector<std::string> test_names(3); > test_names[0] = "Recursive Template"; > test_names[1] = "Template Loop"; > test_names[2] = "Standard pow()"; > > for (unsigned int i=0; i<tests.size(); ++i) > { > clock_t start = clock(); > > tests[i](base_values); > > clock_t stop = clock(); > > const double time_taken = (stop  start)/double(CLOCKS_PER_SEC); > std::cout << test_names[i] << ": " > << std::setprecision(10) > << time_taken > << std::endl; > > } > > return 0; >} > > > > > >The SF.Net email is sponsored by EclipseCon 2004 >Premiere Conference on Open Tools Development and Integration >See the breadth of Eclipse activity. February 35 in Anaheim, CA. >http://www.eclipsecon.org/osdn >_______________________________________________ >Libmeshdevel mailing list >Libmeshdevel@... >https://lists.sourceforge.net/lists/listinfo/libmeshdevel > > 
From: John Peterson <peterson@cf...>  20040209 04:19:36

Benjamin S. Kirk writes: > Help me out here, John... Do you remeber that C++ User's Journal on > scientific programming? I think in near the end of that they describe a > templated pow function that hadles x^int cases intelligently through > repeated squaring... Does that ring a bell? It might be worth having > in utility.h OK, I've done some testing. Apparently pow() is ridiculously slow. I tested two different methods and pow in the attached code. Under O3 optimization, I found the results on a P4 2GHz (in seconds): Raising 50 million double values to the third power. Recursive Template: 0.05 Template Loop: 1.93 Standard pow(): 13.01 Raising to the fifth power, results were: Recursive Template: 0.05 Template Loop: 2.23 Standard pow(): 14.17 Finally, raising to the 10th power (base numbers were small) Recursive Template: 0.06 Template Loop: 2.64 Standard pow(): 13.56 The "Recursive Template" method is the method I posted to the mailing list earlier today. It seems to be the fastest, and as expected, its performance remains about constant no matter how high the exponent becomes. Here is the code (don't forget lm for pow()) I invite you to test it out yourselves... #include <iostream> #include <vector> #include <time.h> #include <iomanip> // Anonymous Namespace to hold useful constants namespace { const unsigned int n_bases = static_cast<unsigned int>(1e7); // Number of tests to perform const unsigned int n_tests = 5; // Raise to an integer power const unsigned int power = 10; }; // Power raising class template <int N> struct Pow { template <int Val> struct In { enum { value = Val * Pow<N1>::template In<Val>::value }; }; }; template <> struct Pow<0> { template <int Val> struct In { enum { value = 1 }; }; }; // Rescursive template function template <int N> double pow(double x) { return x * pow<N1>(x); } template <> double pow<0>(double x) { return 1.; } // Templated Loop function template <int N> double pow2(double x) { const unsigned int l = N / 2; const unsigned int m = N % 2; const double x2 = x*x; double y = 1.0; for (unsigned int i=0; i<l; ++i) { y = y * x2; } return m==0 ? y : y*x; } // I can't figure out how to get the class to use // double as its base value. Therefore, this method // is not feasible now. // void test_class(std::vector<double>& base_values) // { // for (unsigned int n=0; n<n_tests; ++n) // { // for (unsigned int i=0; i<n_bases; ++i) // //std::cout << // Pow<base_values[i]>::In<power> // // << std::endl // ; // } // } void test_recursive_template(std::vector<double>& base_values) { for (unsigned int n=0; n<n_tests; ++n) { for (unsigned int i=0; i<n_bases; ++i) //std::cout << pow<power>(base_values[i]) // << std::endl ; } } void test_template_loop(std::vector<double>& base_values) { for (unsigned int n=0; n<n_tests; ++n) { for (unsigned int i=0; i<n_bases; ++i) //std::cout << pow2<power>(base_values[i]) // << std::endl ; } } void test_standard_pow(std::vector<double>& base_values) { for (unsigned int n=0; n<n_tests; ++n) { for (unsigned int i=0; i<n_bases; ++i) //std::cout << pow(base_values[i], power) // << std::endl ; } } int main() { srand(0); std::vector<double> base_values(n_bases); // Fill the vector with random values between zero and one. for (unsigned int i=0; i<n_bases; ++i) base_values[i] = static_cast<double>(rand()) / static_cast<double>(RAND_MAX); // Convenient typedef for function pointers typedef void (*TestFunction) (std::vector<double>& v); // A vector of function pointers for the functions we will test. std::vector<TestFunction> tests(3); tests[0] = test_recursive_template; tests[1] = test_template_loop; tests[2] = test_standard_pow; // A vector of convenient names std::vector<std::string> test_names(3); test_names[0] = "Recursive Template"; test_names[1] = "Template Loop"; test_names[2] = "Standard pow()"; for (unsigned int i=0; i<tests.size(); ++i) { clock_t start = clock(); tests[i](base_values); clock_t stop = clock(); const double time_taken = (stop  start)/double(CLOCKS_PER_SEC); std::cout << test_names[i] << ": " << std::setprecision(10) << time_taken << std::endl; } return 0; } 
From: John Peterson <peterson@cf...>  20040208 23:28:32

Benjamin S. Kirk writes: > Help me out here, John... Do you remeber that C++ User's Journal on > scientific programming? I think in near the end of that they describe a > templated pow function that hadles x^int cases intelligently through > repeated squaring... Does that ring a bell? It might be worth having > in utility.h It's tough to beat successive squaring to raise to an arbitrary integer power. I would be surprised if pow does not already do that. Those C++ User's Journals are up at school right now, I will take a look through them on Monday. This pair of functions can do it, but I don't if it's faster. Also I did not test it in a library, i.e. do you have to explicitly instantiate every version you want to use? template <int N> double pow(double x) { return x * pow<N1>(x); } template <> double pow<0>(double x) { return 1.; } int main() { std::cout << pow<3>(2.0) << std::endl; return 0; } John 
From: Benjamin S. Kirk <benkirk@cf...>  20040208 22:08:39

I just changed the sqrt calls in the shape function evaluation to constants... As for pow... We might be able to do better than that. In particular, pow is not particularly smart & does not make optimizations where possible. Help me out here, John... Do you remeber that C++ User's Journal on scientific programming? I think in near the end of that they describe a templated pow function that hadles x^int cases intelligently through repeated squaring... Does that ring a bell? It might be worth having in utility.h In the mean time I have added inline pow2, pow3, pow4, pow5 functions in that file. Test it out and let me know if it is fastest (in optimized mode) than the old implementation. Of course, make sure the answers are still the same! Finally, I changed all the tests that look like (a != min(a,b)) to (a>b) I realize that the former is the test I put in the hierarchic implementation... What was I thinking? That code is *really* messy in 3D... I need to look at it again. Ben Daniel Dreyer wrote: > On Sat, 7 Feb 2004, Benjamin S. Kirk wrote: > > >>Thanks for the prompt on this! >> >>Many of the changes in CVS simply changed "Copyright 20022003" with >>"20022004" >> >>However, there are a number of other functional changes as well. I >>guess I can put together an official release this week (comments from >>the devel list?). The biggest job in creating an "official" release is >>in testing the code on all the machines I have access to, and testing >>various configuration options. >> >>I'll try to release 0.4.2 this week if there are no complaints from the >>other developers. >> >>Ben > > > > The SzaboBabuska shape functions are fully functional, and > they can already be used. However, they are not yet _truly_ > good: only the 1D shapes are coded in Horner scheme. > > The 2D shapes, honestly, aren't perfect (yet): many > calls to sqrt(), pow(x,n) with n even 5 or so... > That means: the 2D shapes (namely only on TRI6) may be prone > to cancellation issues, and very likely slow! > > I'd suggest to: > >  either change everything _prior_ to releasing 0.4.2, > namely all the sqrt to evaluated floats, and the > pow's to hornerlike. Since Steffen/Hendrik did these shapes, > i'd suggest they work on that. They have to say how long > this may take. > >  or postpone this to 0.4.3, and put some comment in the release text > that says something like "SZABAB is functional (compiles fine), however > not really finished". > > Either way, i suggest that the pows and sqrts should go sooner > or later, either before or after 0.4.2. > > >  Sorry for being that pedantic, but since libmesh is > open source, we better have an eye on the quality of > the code. > > Daniel 
From: Benjamin S. Kirk <benkirk@cf...>  20040208 20:35:37

Some time ago John said that a compiler configuration script would be useful when building test programs... The goal here is to not require a Makefile to build simple code. I agreed, but until now didn't have much use for such a thing, so I did nothing... :( I just added libmeshconfig. (Actually, contrib/bin/libmeshconfig.in, which is processed by configure...) It is pretty simple, please improve it in any way you see fit... Here is the idea: c++ o foo foo.C `libmeshconfig cxxflags include ldflags` will build foo. The big deal is that all the paths & flags are properly set, without the intervention of a Makefile. If only PETSc came with one of these... Ben 
From: Benjamin S. Kirk <benkirk@cf...>  20040208 00:54:21

Thanks for the prompt on this! Many of the changes in CVS simply changed "Copyright 20022003" with "20022004" However, there are a number of other functional changes as well. I guess I can put together an official release this week (comments from the devel list?). The biggest job in creating an "official" release is in testing the code on all the machines I have access to, and testing various configuration options. I'll try to release 0.4.2 this week if there are no complaints from the other developers. Ben Denis Danilov wrote: >Hi, > >According CVS i see many things have been changed in >comparison with libmesh0.4.1. Is there any plans about >next stable release of libmesh in near future? > >Best regards, >Denis > > > >The SF.Net email is sponsored by EclipseCon 2004 >Premiere Conference on Open Tools Development and Integration >See the breadth of Eclipse activity. February 35 in Anaheim, CA. >http://www.eclipsecon.org/osdn >_______________________________________________ >Libmeshusers mailing list >Libmeshusers@... >https://lists.sourceforge.net/lists/listinfo/libmeshusers > > 
From: John Peterson <peterson@cf...>  20040205 23:11:10

Daniel Dreyer writes: > suggested by John. But why distinguish between C0_FE and C1_FE, > then where should go the pFEM like SZABAB go in? Why not this: > > FE < InfFE > FE < ScalarFE > FE < VectorFE > > and ScalarFE is just templated as before FE<d,family>? Actually I was thinking about this so I am glad somebody else brought it up. This hierarchy seems to be a more general distinction. Subclasses of the ScalarFE class are free to be C0, C1, discontinuous, etc. John 
From: KIRK, BENJAMIN (JSCEG) (NASA) <benjamin.kirk1@na...>  20040203 14:29:48

I don't really see any problem with that. I think the implementation could be handled just like everything else in the FE class, i.e. through template specialization. The current (scalar) Finite elements could continue to work as they are, they would just operate on the subset phi[0][i][qp] and get_phi() would return a reference to phi[0]. I think this could easily handle the RaviartThomas part, and they should be easy to implement. I think the Nedelec part might be more complicated, but it should also fit into this framework. On a related note, maybe the best way to handle this is to expand the FEBase/FE setup. Perhaps a pure virtual base class FE with a range of classes that inherit from it: FE<C0_FE<Lagrange... FE<InfFE... FE<C1_FE<Hermite... FE<DFE<Discontinuous... (this is where a whiteboard is handy). So, we might reexamine the FE structure and think about improving it. There is a fair amount of code in there, though... Thoughts? Ben Original Message From: John Peterson [mailto:peterson@...] Sent: Sunday, February 01, 2004 9:36 PM To: KIRK, BENJAMIN (JSCEG) (NASA) Subject: vectorvalued shape functions I've recently been asked about programming the RaviartThomasNedelec elements in libMesh by someone I met at the Manchester conference from this summer. Although I don't have the reference yet (J. C. Nedelec, 1980, "Mixed finite elements in R3." Numerical Math. 35:315341) I gather that they are vectorvalued shape functions, i.e. when you dot the governing PDE with a vector shape function v, v is no longer simply (phi_i,0,0), (0,phi_i,0), (0,0,phi_i), the components are dependent to enforce some property such as divergence free, etc. Wolfgang has a discussion on matrix assembly for these types of elements in deal.ii here: http://gaia.iwr.uniheidelberg.de/~deal/developer/documentation.html (search for assembling matrices on the page) it seems like a very generic idea since the scalarvalued shape functions now are simply a special case where only one component of the vector is nonzero. We might add a function called "is_primitive" which returns whether or not a particular FE instantiation has vectorvalued shape functions. Or we could do something fancier with virtual functions... The functions get_phi, get_dphi would fail for nonprimitive elements and continue working for the current elements. For nonprimitive elements you would need to call for example: get_phi_vector() which would return a reference to a vector of vectors of vectors. Then use #defines as normal: #define phix(i) phi_values[0][i][qp] #define phiy(i) phi_values[1][i][qp] #define phiz(i) phi_values[2][i][qp] Do you see any big problems with that? It might lead to more users for the library since many people solving Maxwell's equation often use H(curl) conforming elements which fit into this framework. John 