From: Andrew E S. <and...@gm...> - 2012-06-27 18:25:58
|
If I understand libmesh correctly, if you are solving a system that has a vector solution (i.e., velocity) then you need to add variables such as "vx" and "vy" for each component. Then solve the system using submatrices and vectors etc. Therefore, the libmesh::System add_vector command is only for storing supplement data. Is this correct libmesh logic? If I am correct above, then how do a write an initialization function for such a problem? I attached a sample program below that doesn't work because the initial_velocity function is called for both and the first value of the output vector is assigned to both vx and vy. How do I assign different values for the different directions? Thanks, Andrew // Standard library includes #include <iostream> #include <algorithm> #include <math.h> #include <stdio.h> #include <string> using std::string; // libMesh includes #include "libmesh.h" #include "mesh.h" #include "mesh_generation.h" #include "equation_systems.h" #include "nonlinear_implicit_system.h" #include "transient_system.h" #include "explicit_system.h" #include "analytic_function.h" #include "exodusII_io.h" // Bring in everything from the libMesh namespace using namespace libMesh; // A function for initializing, this d void initial_velocity(DenseVector<Number>& output, const Point&, const Real){ output(0) = 1; // x-dir // output(1) = 2; // y-dir (this causes an error, but this is what seems logical) } void init_cd (EquationSystems& es, const std::string&){ TransientNonlinearImplicitSystem& system = es.get_system<TransientNonlinearImplicitSystem>("test"); AnalyticFunction<Number> fobj(initial_velocity); system.project_solution(&fobj); } // The main function int main (int argc, char** argv){ // Initialize libraries LibMeshInit init (argc, argv); // Generate a mesh Mesh mesh; MeshTools::Generation::build_square(mesh, 1, 1, -1., 1, -1, 1, QUAD4); mesh.all_first_order(); // Create an equation system EquationSystems eq_sys(mesh); // Create a system and add two variables (x and y velocity components) eq_sys.add_system<TransientNonlinearImplicitSystem>("test"); eq_sys.get_system("test").add_variable("vx", FIRST); eq_sys.get_system("test").add_variable("vy", FIRST); // Attach initilization function eq_sys.get_system("test").attach_init_function(init_cd); // Initilize the equation system storing the data eq_sys.get_system("test").init(); // Export the solution ExodusII_IO(mesh).write_equation_systems("example3b.ex2", eq_sys); // Done return 0; } -- Andrew E. Slaughter, PhD and...@gm... Materials Process Design and Control Laboratory Sibley School of Mechanical and Aerospace Engineering 169 Frank H. T. Rhodes Hall Cornell University Ithaca, NY 14853-3801 (607) 229-1829 http://aeslaughter.github.com/ |
From: John P. <jwp...@gm...> - 2012-06-27 18:30:37
|
On Wed, Jun 27, 2012 at 12:25 PM, Andrew E Slaughter <and...@gm...> wrote: > > // A function for initializing, this d > void initial_velocity(DenseVector<Number>& output, const Point&, const > Real){ > output(0) = 1; // x-dir > // output(1) = 2; // y-dir (this causes an error, but this is what seems > logical) > } I think you are on the right track here. It's just that output is not resized for you, so output.resize(2); and try again? -- John |
From: Andrew E S. <and...@gm...> - 2012-06-27 18:44:58
|
I was playing around with that, but the function still gets called for each variable. Essentially, doing the calculations twice or three times in 3D. This seems far to inefficient. Is there no way around this? void initial_velocity(DenseVector<Number>& output, const Point&, const Real){ output.resize(2); output(0) = 1; // x-direction velocity output(1) = 2; // y-direction velocity printf("hello\n"); // you get 8 of these (4 is expected) } On Wed, Jun 27, 2012 at 2:30 PM, John Peterson <jwp...@gm...> wrote: > On Wed, Jun 27, 2012 at 12:25 PM, Andrew E Slaughter > <and...@gm...> wrote: > > > > // A function for initializing, this d > > void initial_velocity(DenseVector<Number>& output, const Point&, const > > Real){ > > output(0) = 1; // x-dir > > // output(1) = 2; // y-dir (this causes an error, but this is what seems > > logical) > > } > > I think you are on the right track here. It's just that output is not > resized for you, so output.resize(2); and try again? > > -- > John > -- Andrew E. Slaughter, PhD and...@gm... Materials Process Design and Control Laboratory Sibley School of Mechanical and Aerospace Engineering 169 Frank H. T. Rhodes Hall Cornell University Ithaca, NY 14853-3801 (607) 229-1829 http://aeslaughter.github.com/ |
From: Paul T. B. <ptb...@gm...> - 2012-06-27 18:55:13
|
On Wed, Jun 27, 2012 at 1:44 PM, Andrew E Slaughter < and...@gm...> wrote: > I was playing around with that, but the function still gets called for each > variable. Essentially, doing the calculations twice or three times in 3D. > This seems far to inefficient. Is there no way around this? > Internally, libMesh is calling the component method of AnalyticFunction. It is virtual, so you can overload it to not compute the whole vector each time. HTH, Paul |
From: John P. <jwp...@gm...> - 2012-06-27 19:01:01
|
On Wed, Jun 27, 2012 at 12:55 PM, Paul T. Bauman <ptb...@gm...> wrote: > > > On Wed, Jun 27, 2012 at 1:44 PM, Andrew E Slaughter > <and...@gm...> wrote: >> >> I was playing around with that, but the function still gets called for >> each >> variable. Essentially, doing the calculations twice or three times in 3D. >> This seems far to inefficient. Is there no way around this? > > > Internally, libMesh is calling the component method of AnalyticFunction. It > is virtual, so you can overload it to not compute the whole vector each > time. Paul's right. Here's the comment for component() from function_base.h * Subclasses are recommended to overload this, since the default * implementation is based on a vector evaluation, which is usually * unnecessarily inefficient. I guess if your initial condition is super simple (like return 1;) you can avoid implementing component() and just do extra work. On the other hand, you can implement component() with a big if-block if computing the initial condition is very expensive... if (i==0) return 1; else return 2; -- John |
From: Andrew E S. <and...@gm...> - 2012-06-27 20:12:07
|
If I understand correctly, if my initialization is expensive then it should be coded directly into the component member of a class derived from AnalyticFunction, as follows. Then simply use the project_solution and skip the add_initial_function(...) and init() commands for the system. Is this correct? class AnalyticFunctionTest : public AnalyticFunction<Number>{ public: AnalyticFunctionTest(){ _initialized = true; } Number component(unsigned int i, const Point& p, Real time = 0){ if (i == 0){ return 1; } else { return 2; } } }; I tried using this, but I get an error b/c I believe that AnalyticFunction needs a default constructor ( http://www.cplusplus.com/doc/tutorial/inheritance/), but it doesn't have one. example3b.cpp:41:25: error: no matching function for call to ‘libMesh::AnalyticFunction<double>::AnalyticFunction()’ So, then how do you do what you propose? On Wed, Jun 27, 2012 at 3:00 PM, John Peterson <jwp...@gm...> wrote: > On Wed, Jun 27, 2012 at 12:55 PM, Paul T. Bauman <ptb...@gm...> > wrote: > > > > > > On Wed, Jun 27, 2012 at 1:44 PM, Andrew E Slaughter > > <and...@gm...> wrote: > >> > >> I was playing around with that, but the function still gets called for > >> each > >> variable. Essentially, doing the calculations twice or three times in > 3D. > >> This seems far to inefficient. Is there no way around this? > > > > > > Internally, libMesh is calling the component method of AnalyticFunction. > It > > is virtual, so you can overload it to not compute the whole vector each > > time. > > Paul's right. Here's the comment for component() from function_base.h > > * Subclasses are recommended to overload this, since the default > * implementation is based on a vector evaluation, which is usually > * unnecessarily inefficient. > > I guess if your initial condition is super simple (like return 1;) you > can avoid implementing component() and just do extra work. > > On the other hand, you can implement component() with a big if-block > if computing the initial condition is very expensive... > > if (i==0) > return 1; > else > return 2; > > -- > John > -- Andrew E. Slaughter, PhD and...@gm... Materials Process Design and Control Laboratory Sibley School of Mechanical and Aerospace Engineering 169 Frank H. T. Rhodes Hall Cornell University Ithaca, NY 14853-3801 (607) 229-1829 http://aeslaughter.github.com/ |
From: John P. <jwp...@gm...> - 2012-06-27 20:25:21
|
On Wed, Jun 27, 2012 at 2:11 PM, Andrew E Slaughter <and...@gm...> wrote: > If I understand correctly, if my initialization is expensive then it should > be coded directly into the component member of a class derived from > AnalyticFunction, as follows. Then simply use the project_solution and skip > the add_initial_function(...) and init() commands for the system. Is this > correct? > > class AnalyticFunctionTest : public AnalyticFunction<Number>{ > public: > AnalyticFunctionTest(){ > _initialized = true; > } > > Number component(unsigned int i, const Point& p, Real time = 0){ > if (i == 0){ > return 1; > } else { > return 2; > } > } > }; > > I tried using this, but I get an error b/c I believe that AnalyticFunction > needs a default constructor > (http://www.cplusplus.com/doc/tutorial/inheritance/), but it doesn't have > one. > > example3b.cpp:41:25: error: no matching function for call to > ‘libMesh::AnalyticFunction<double>::AnalyticFunction()’ > > So, then how do you do what you propose? I think you want to write a forwarding constructor for your class that takes a pointer-to-function and constructs the parent by passing it along: class AnalyticFunctionTest : public AnalyticFunction<Number>{ public: AnalyticFunctionTest(void fptr(DenseVector<Output>& output, const Point& p, const Real time)) : AnalyticFunction(fptr) // <-- construct the parent { _initialized = true; } -- John |
From: Andrew E S. <and...@gm...> - 2012-06-27 20:30:23
|
But, if I am programming the component member to perform the initialization, I don't need a function pointer, so passing one in is meaningless. Is this just an unfortunate consequence or is the pointer actually going to be used in some other capacity? Thanks, Andrew On Wed, Jun 27, 2012 at 4:24 PM, John Peterson <jwp...@gm...> wrote: > On Wed, Jun 27, 2012 at 2:11 PM, Andrew E Slaughter > <and...@gm...> wrote: > > If I understand correctly, if my initialization is expensive then it > should > > be coded directly into the component member of a class derived from > > AnalyticFunction, as follows. Then simply use the project_solution and > skip > > the add_initial_function(...) and init() commands for the system. Is this > > correct? > > > > class AnalyticFunctionTest : public AnalyticFunction<Number>{ > > public: > > AnalyticFunctionTest(){ > > _initialized = true; > > } > > > > Number component(unsigned int i, const Point& p, Real time = 0){ > > if (i == 0){ > > return 1; > > } else { > > return 2; > > } > > } > > }; > > > > I tried using this, but I get an error b/c I believe that > AnalyticFunction > > needs a default constructor > > (http://www.cplusplus.com/doc/tutorial/inheritance/), but it doesn't > have > > one. > > > > example3b.cpp:41:25: error: no matching function for call to > > ‘libMesh::AnalyticFunction<double>::AnalyticFunction()’ > > > > So, then how do you do what you propose? > > I think you want to write a forwarding constructor for your class that > takes a pointer-to-function and constructs the parent by passing it > along: > > class AnalyticFunctionTest : public AnalyticFunction<Number>{ > public: > AnalyticFunctionTest(void fptr(DenseVector<Output>& output, > const Point& p, > const Real time)) > : AnalyticFunction(fptr) // <-- construct the parent > { > _initialized = true; > } > > > -- > John > -- Andrew E. Slaughter, PhD and...@gm... Materials Process Design and Control Laboratory Sibley School of Mechanical and Aerospace Engineering 169 Frank H. T. Rhodes Hall Cornell University Ithaca, NY 14853-3801 (607) 229-1829 http://aeslaughter.github.com/ |
From: John P. <jwp...@gm...> - 2012-06-27 20:41:18
|
On Wed, Jun 27, 2012 at 2:30 PM, Andrew E Slaughter <and...@gm...> wrote: > But, if I am programming the component member to perform the initialization, > I don't need a function pointer, so passing one in is meaningless. Is this > just an unfortunate consequence or is the pointer actually going to be used > in some other capacity? Hmm... you're right. I think if you plan to overload the component() it may be better to inherit from AnalyticFunction's base class (FunctionBase) instead. If you do that, you'll have to provide operator()(Point, Real) operator()(Point, Real, DenseVector) clone() and component() The first 3 are pure virtual so you have to overload them, but I'm not sure specifically which of the operator() functions are actually called (hopefully only component() is called?) Paul, do you have an example code where you actually overload component() since you suggested it? -- John |
From: Paul T. B. <ptb...@gm...> - 2012-06-27 21:52:59
|
On Wed, Jun 27, 2012 at 3:40 PM, John Peterson <jwp...@gm...> wrote: > > Paul, do you have an example code where you actually overload > component() since you suggested it? > Yeah, I would suggest deriving from FunctionBase directly. Check out boundary_function.h in examples/vector_fe/vector_fe_ex2 (this will be in trunk within the past couple of days) for an example of this. AFAICT, libMesh internally only uses component, but I could be wrong. Certiainly, this is what is called in the projection routines. HTH, Paul |
From: Andrew E S. <and...@gm...> - 2012-06-27 23:36:05
|
Paul, I just checkout out the trunk with the following and the folder you specify does not exist, only examples/vector_fe/vector_fe_ex1.C svn checkout https://libmesh.svn.sourceforge.net/svnroot/libmesh/trunk/libmesh On Wed, Jun 27, 2012 at 5:52 PM, Paul T. Bauman <ptb...@gm...> wrote: > On Wed, Jun 27, 2012 at 3:40 PM, John Peterson <jwp...@gm...> > wrote: > > > > > Paul, do you have an example code where you actually overload > > component() since you suggested it? > > > > Yeah, I would suggest deriving from FunctionBase directly. Check out > boundary_function.h in examples/vector_fe/vector_fe_ex2 (this will be in > trunk within the past couple of days) for an example of this. AFAICT, > libMesh internally only uses component, but I could be wrong. Certiainly, > this is what is called in the projection routines. > > HTH, > > Paul > > ------------------------------------------------------------------------------ > Live Security Virtual Conference > Exclusive live event will cover all the ways today's security and > threat landscape has changed and how IT managers can respond. Discussions > will include endpoint security, mobile security and the latest in malware > threats. http://www.accelacomm.com/jaw/sfrnl04242012/114/50122263/ > _______________________________________________ > Libmesh-users mailing list > Lib...@li... > https://lists.sourceforge.net/lists/listinfo/libmesh-users > -- Andrew E. Slaughter, PhD and...@gm... Materials Process Design and Control Laboratory Sibley School of Mechanical and Aerospace Engineering 169 Frank H. T. Rhodes Hall Cornell University Ithaca, NY 14853-3801 (607) 229-1829 http://aeslaughter.github.com/ |
From: Paul T. B. <ptb...@gm...> - 2012-06-28 00:31:10
|
On Wed, Jun 27, 2012 at 6:35 PM, Andrew E Slaughter < and...@gm...> wrote: > I just checkout out the trunk with the following and the folder you > specify does not exist, only examples/vector_fe/vector_fe_ex1.C > *facepalm* Sorry, that patch is still being discussed on libmesh-devel. You can grab the function here (I also include the function that it calls for posterity): http://users.ices.utexas.edu/~pbauman/bound_func/ Sorry about that. Paul |
From: Andrew E S. <and...@gm...> - 2012-06-28 01:20:50
|
Paul, Thanks for linking the function, that is where I will head with my initialization. I will post my implementation when I am done in case anyone needs this in the future. One question, where is libmesh_not_implemented() defined, I can't seem to find any documentation on it, but it looks useful? Thanks, Andrew On Wed, Jun 27, 2012 at 8:31 PM, Paul T. Bauman <ptb...@gm...> wrote: > On Wed, Jun 27, 2012 at 6:35 PM, Andrew E Slaughter < > and...@gm...> wrote: > >> I just checkout out the trunk with the following and the folder you >> specify does not exist, only examples/vector_fe/vector_fe_ex1.C >> > > *facepalm* > > Sorry, that patch is still being discussed on libmesh-devel. You can grab > the function here (I also include the function that it calls for > posterity): http://users.ices.utexas.edu/~pbauman/bound_func/ > > Sorry about that. > > Paul > -- Andrew E. Slaughter, PhD and...@gm... Materials Process Design and Control Laboratory Sibley School of Mechanical and Aerospace Engineering 169 Frank H. T. Rhodes Hall Cornell University Ithaca, NY 14853-3801 (607) 229-1829 http://aeslaughter.github.com/ |
From: Andrew E S. <and...@gm...> - 2012-06-28 01:40:41
|
Here is my test problem with an initialization class that inherits from FunctionBase. Thanks for the help! -Andrew // Standard library includes #include <iostream> #include <algorithm> #include <math.h> #include <stdio.h> #include <string> using std::string; // libMesh includes #include "libmesh.h" #include "mesh.h" #include "mesh_generation.h" #include "equation_systems.h" #include "nonlinear_implicit_system.h" #include "transient_system.h" #include "explicit_system.h" #include "analytic_function.h" #include "exodusII_io.h" // Bring in everything from the libMesh namespace using namespace libMesh; // A class for initilizing velocity template <typename Output = Number> class EqInit : public FunctionBase<Output>{ public: // Returns the velocity by its components virtual Number component(unsigned int index, const Point& p, Real time = 0){ if (index == 0){ return 1; // add some intensive calculation here } else { return 2; // add some intensive calculation here } } // Do not use the scalar output version of the operator() Output operator() (const Point& p, const Real t = 0.){ libmesh_not_implemented(); } // Use the component member to build the vector output version of operator() void operator() (const Point& p, const Real t, DenseVector<Output>& output){ output.resize(2); output(0) = component(0,p,t); output(1) = component(1,p,t); return; } // Create a copy of the class virtual AutoPtr<FunctionBase<Output> > clone () const{ return AutoPtr<FunctionBase<Output> > (new EqInit<Output>()); } }; // The main function int main (int argc, char** argv){ // Initialize libraries LibMeshInit init (argc, argv); // Generate a mesh Mesh mesh; MeshTools::Generation::build_square(mesh, 1, 1, -1., 1, -1, 1, QUAD4); mesh.all_first_order(); // Create an equation system EquationSystems eq_sys(mesh); // Create a system and add two variables (x and y velocity components) eq_sys.add_system<TransientNonlinearImplicitSystem>("test"); eq_sys.get_system("test").add_variable("vx", FIRST); eq_sys.get_system("test").add_variable("vy", FIRST); // Initilize the equation system eq_sys.get_system("test").init(); // Project the velocity using the EqInit class EqInit<Number> init_velocity; eq_sys.get_system("test").project_solution(&init_velocity); // Test that values are correct for each point on the only element DenseVector<Number> v(2); Elem* elem = mesh.elem(0); for (unsigned int i = 0; i < elem->n_nodes(); i++){ Point& p = elem->point(i); v(0) = eq_sys.get_system("test").point_value(0,p); v(1) = eq_sys.get_system("test").point_value(1,p); printf("Point %g,%g; v = %f,%f\n", p(0),p(1),v(0),v(1)); } // Done return 0; } On Wed, Jun 27, 2012 at 9:20 PM, Andrew E Slaughter < and...@gm...> wrote: > Paul, > Thanks for linking the function, that is where I will head with > my initialization. I will post my implementation when I am done in case > anyone needs this in the future. > > One question, where is libmesh_not_implemented() defined, I can't seem to > find any documentation on it, but it looks useful? > > Thanks, > Andrew > > > On Wed, Jun 27, 2012 at 8:31 PM, Paul T. Bauman <ptb...@gm...>wrote: > >> On Wed, Jun 27, 2012 at 6:35 PM, Andrew E Slaughter < >> and...@gm...> wrote: >> >>> I just checkout out the trunk with the following and the folder you >>> specify does not exist, only examples/vector_fe/vector_fe_ex1.C >>> >> >> *facepalm* >> >> Sorry, that patch is still being discussed on libmesh-devel. You can grab >> the function here (I also include the function that it calls for >> posterity): http://users.ices.utexas.edu/~pbauman/bound_func/ >> >> Sorry about that. >> >> Paul >> > > > > -- > Andrew E. Slaughter, PhD > and...@gm... > > Materials Process Design and Control Laboratory > Sibley School of Mechanical and Aerospace Engineering > 169 Frank H. T. Rhodes Hall > Cornell University > Ithaca, NY 14853-3801 > (607) 229-1829 > http://aeslaughter.github.com/ > > -- Andrew E. Slaughter, PhD and...@gm... Materials Process Design and Control Laboratory Sibley School of Mechanical and Aerospace Engineering 169 Frank H. T. Rhodes Hall Cornell University Ithaca, NY 14853-3801 (607) 229-1829 http://aeslaughter.github.com/ |
From: Paul T. B. <ptb...@gm...> - 2012-06-28 01:44:18
|
On Wed, Jun 27, 2012 at 8:20 PM, Andrew E Slaughter < and...@gm...> wrote: > Paul, > Thanks for linking the function, that is where I will head with > my initialization. I will post my implementation when I am done in case > anyone needs this in the future. > No problem, glad it was helpful. > One question, where is libmesh_not_implemented() defined, I can't seem to > find any documentation on it, but it looks useful? > libmesh_common.h You'll find many other useful macros there as well. Best, Paul |