From: Manav B. <bha...@gm...> - 2013-09-20 18:34:26
|
Hi, Have there been any changes to System::project_solution() with FEMFunctionBase over the last 2-3 months? I had a piece of code that was working fine with it, but now is throwing an error at line 2339 in sytem_projection.C: const unsigned int n_qp = qsiderule.n_points(); It seems like the FEMContext object (context) initialized at line 1843 in system_projection.C sets up its element_qrule and side_qrule with zero points, but edge_qrule has 2 points. Perhaps an initialization needs to be added? I would appreciate any help. Thanks, Manav |
From: Paul T. B. <ptb...@gm...> - 2013-09-20 18:45:45
|
On Fri, Sep 20, 2013 at 1:34 PM, Manav Bhatia <bha...@gm...> wrote: > > Have there been any changes to System::project_solution() with > FEMFunctionBase over the last 2-3 months? > Current master shows FEMFunctionBase was touched 5 months ago. > I had a piece of code that was working fine with it, but now is throwing > an error at line 2339 in sytem_projection.C: > Is this with the current git master? 2339 in system_projection.C is not a libmesh_error()... What's the error you're getting? > const unsigned int n_qp = qsiderule.n_points(); > > It seems like the FEMContext object (context) initialized at line 1843 in > system_projection.C sets up its element_qrule and side_qrule with zero > points, but edge_qrule has 2 points. > This makes me a little nervous because we pushed API updates to FEMContext, but it shouldn't have affected any of that (unless I screwed something up obviously). Again, version of libMesh and the error your getting? stacktrace? Thanks, Paul |
From: Manav B. <bha...@gm...> - 2013-09-20 19:05:19
|
On Fri, Sep 20, 2013 at 2:45 PM, Paul T. Bauman <ptb...@gm...> wrote: > On Fri, Sep 20, 2013 at 1:34 PM, Manav Bhatia <bha...@gm...>wrote: > >> >> Have there been any changes to System::project_solution() with >> FEMFunctionBase over the last 2-3 months? >> > > Current master shows FEMFunctionBase was touched 5 months ago. > > >> I had a piece of code that was working fine with it, but now is >> throwing >> an error at line 2339 in sytem_projection.C: >> > > Is this with the current git master? 2339 in system_projection.C is not a > libmesh_error()... What's the error you're getting? > Sorry I was not too specific. I am on the git master, and line 2339 is a call to qsiderule.n_points(). It has an assertion of libmesh_assert(!_points.empyt()), where the code throws an exception. The stacktrace is (my internet machine is different from my coding machine, so I am typing this out): 0: 0 libMesh::print_trace 1: 1 libMesh::MacroFunctions::report_error 2: 2 QBase::n_points 3: 3 libMesh::ProjectFEMSolution::operator() 4: 4 libMesh::Threads::parallel_for<>... 5: 5 libMesh::System::project_vector 6: 6 libMesh::System::project_solution 7: 7 call from application code main Please let me know if you would need more information. > const unsigned int n_qp = qsiderule.n_points(); >> >> It seems like the FEMContext object (context) initialized at line 1843 >> in >> system_projection.C sets up its element_qrule and side_qrule with zero >> points, but edge_qrule has 2 points. >> > > This makes me a little nervous because we pushed API updates to > FEMContext, but it shouldn't have affected any of that (unless I screwed > something up obviously). > > Again, version of libMesh and the error your getting? stacktrace? > > Thanks, > > Paul > > |
From: Roy S. <roy...@ic...> - 2013-09-20 20:07:35
|
On Fri, 20 Sep 2013, Manav Bhatia wrote: > Sorry I was not too specific. I am on the git master, and line 2339 is a > call to qsiderule.n_points(). It has an assertion of > libmesh_assert(!_points.empyt()), where the code throws an exception. I think Paul's (and my) confusion was a simple typo - you're referring to line 2239, not 2339. --- Roy |
From: Manav B. <bha...@gm...> - 2013-09-20 20:09:28
|
Yes, it is 2239. I am sorry about the typo and the resulting confusion. -Manav On Fri, Sep 20, 2013 at 4:07 PM, Roy Stogner <roy...@ic...>wrote: > > On Fri, 20 Sep 2013, Manav Bhatia wrote: > > Sorry I was not too specific. I am on the git master, and line 2339 is a >> call to qsiderule.n_points(). It has an assertion of >> libmesh_assert(!_points.empyt(**)), where the code throws an exception. >> > > I think Paul's (and my) confusion was a simple typo - you're referring > to line 2239, not 2339. > --- > Roy > |
From: Manav B. <bha...@gm...> - 2013-09-20 20:21:16
|
I was comparing the code with AssemblyContributions in fem_system.C, and see that context.side_fe_reinit() has not been called before the call to get qsiderule.n_points(). Should that be added to the for loop on element sides at line 2248 and move the get_side_qrule method call (line 2238) after this call? -Manav On Fri, Sep 20, 2013 at 4:09 PM, Manav Bhatia <bha...@gm...> wrote: > Yes, it is 2239. I am sorry about the typo and the resulting confusion. > -Manav > > > On Fri, Sep 20, 2013 at 4:07 PM, Roy Stogner <roy...@ic...>wrote: > >> >> On Fri, 20 Sep 2013, Manav Bhatia wrote: >> >> Sorry I was not too specific. I am on the git master, and line 2339 is a >>> call to qsiderule.n_points(). It has an assertion of >>> libmesh_assert(!_points.empyt(**)), where the code throws an exception. >>> >> >> I think Paul's (and my) confusion was a simple typo - you're referring >> to line 2239, not 2339. >> --- >> Roy >> > > |
From: Manav B. <bha...@gm...> - 2013-09-20 20:24:31
|
Just found the the call to side_fe_reinit() in this for loop. It seems like the entire block from lines 2238 to 2247 should be moved to after the context.side_fe_reinit() ? -Manav On Fri, Sep 20, 2013 at 4:21 PM, Manav Bhatia <bha...@gm...> wrote: > I was comparing the code with AssemblyContributions in fem_system.C, and > see that context.side_fe_reinit() has not been called before the call to > get qsiderule.n_points(). Should that be added to the for loop on element > sides at line 2248 and move the get_side_qrule method call (line > 2238) after this call? > > -Manav > > On Fri, Sep 20, 2013 at 4:09 PM, Manav Bhatia <bha...@gm...>wrote: > >> Yes, it is 2239. I am sorry about the typo and the resulting confusion. >> -Manav >> >> >> On Fri, Sep 20, 2013 at 4:07 PM, Roy Stogner <roy...@ic...>wrote: >> >>> >>> On Fri, 20 Sep 2013, Manav Bhatia wrote: >>> >>> Sorry I was not too specific. I am on the git master, and line 2339 is a >>>> call to qsiderule.n_points(). It has an assertion of >>>> libmesh_assert(!_points.empyt(**)), where the code throws an exception. >>>> >>> >>> I think Paul's (and my) confusion was a simple typo - you're referring >>> to line 2239, not 2339. >>> --- >>> Roy >>> >> >> > |
From: Paul T. B. <ptb...@gm...> - 2013-09-20 22:21:58
|
Hmm. I can't reproduce this. I have a code that uses FEMFunction for projections and I bumped up the basis order and cannot trip the error you are. In fact, I manually verified that n_points() returns the correct value on that line. Can you make a test case that trips the error? Are you using a LAGRANGE basis or some other one? What order? On Fri, Sep 20, 2013 at 3:24 PM, Manav Bhatia <bha...@gm...> wrote: > Just found the the call to side_fe_reinit() in this for loop. It seems > like the entire block from lines 2238 to 2247 should be moved to after the > context.side_fe_reinit() ? > > -Manav > On Fri, Sep 20, 2013 at 4:21 PM, Manav Bhatia <bha...@gm...>wrote: > >> I was comparing the code with AssemblyContributions in fem_system.C, and >> see that context.side_fe_reinit() has not been called before the call to >> get qsiderule.n_points(). Should that be added to the for loop on element >> sides at line 2248 and move the get_side_qrule method call (line >> 2238) after this call? >> >> -Manav >> >> On Fri, Sep 20, 2013 at 4:09 PM, Manav Bhatia <bha...@gm...>wrote: >> >>> Yes, it is 2239. I am sorry about the typo and the resulting confusion. >>> -Manav >>> >>> >>> On Fri, Sep 20, 2013 at 4:07 PM, Roy Stogner <roy...@ic...>wrote: >>> >>>> >>>> On Fri, 20 Sep 2013, Manav Bhatia wrote: >>>> >>>> Sorry I was not too specific. I am on the git master, and line 2339 is a >>>>> call to qsiderule.n_points(). It has an assertion of >>>>> libmesh_assert(!_points.empyt(**)), where the code throws an >>>>> exception. >>>>> >>>> >>>> I think Paul's (and my) confusion was a simple typo - you're referring >>>> to line 2239, not 2339. >>>> --- >>>> Roy >>>> >>> >>> >> > |
From: Manav B. <bha...@gm...> - 2013-09-21 03:23:03
|
Hi Paul, Did you try this in 2D or 3D? I do not get this error in 2D, but it shows up in 3D. The following code is able to reproduce the issue at my end. int main (int argc, char* const argv[]) { // Initialize libMesh. LibMeshInit init (argc, argv); { ParallelMesh m(init.comm()); MeshTools::Generation::build_cube(m, 5,5,5, 0., 1., 0., 1., 0., 1., HEX8); EquationSystems es(m); ExplicitSystem& s = es.add_system<ExplicitSystem>("mys"); FEType fe; s.add_variable("v", fe); es.init(); class Func: public FEMFunctionBase<Real> { public: virtual AutoPtr<FEMFunctionBase<Real> > clone () const { return AutoPtr<FEMFunctionBase<Real> > (new Func); }; virtual Real operator() (const FEMContext&, const Point& p, const Real time = 0.) { return 0.; }; virtual void operator() (const FEMContext&, const Point& p, const Real time, DenseVector<Real>& output) { output.zero(); }; }; Func ff; s.project_solution(&ff); } return 0; } Please let me know if you see the same behavior. Maybe I am doing something wrong? Thanks, Manav On Sep 20, 2013, at 6:21 PM, Paul T. Bauman <ptb...@gm...> wrote: > Hmm. I can't reproduce this. I have a code that uses FEMFunction for projections and I bumped up the basis order and cannot trip the error you are. In fact, I manually verified that n_points() returns the correct value on that line. Can you make a test case that trips the error? Are you using a LAGRANGE basis or some other one? What order? > > > On Fri, Sep 20, 2013 at 3:24 PM, Manav Bhatia <bha...@gm...> wrote: > Just found the the call to side_fe_reinit() in this for loop. It seems like the entire block from lines 2238 to 2247 should be moved to after the context.side_fe_reinit() ? > > -Manav > On Fri, Sep 20, 2013 at 4:21 PM, Manav Bhatia <bha...@gm...> wrote: > I was comparing the code with AssemblyContributions in fem_system.C, and see that context.side_fe_reinit() has not been called before the call to get qsiderule.n_points(). Should that be added to the for loop on element sides at line 2248 and move the get_side_qrule method call (line 2238) after this call? > > -Manav > > On Fri, Sep 20, 2013 at 4:09 PM, Manav Bhatia <bha...@gm...> wrote: > Yes, it is 2239. I am sorry about the typo and the resulting confusion. > -Manav > > > On Fri, Sep 20, 2013 at 4:07 PM, Roy Stogner <roy...@ic...> wrote: > > On Fri, 20 Sep 2013, Manav Bhatia wrote: > > Sorry I was not too specific. I am on the git master, and line 2339 is a > call to qsiderule.n_points(). It has an assertion of > libmesh_assert(!_points.empyt()), where the code throws an exception. > > I think Paul's (and my) confusion was a simple typo - you're referring > to line 2239, not 2339. > --- > Roy > > > > |
From: Paul T. B. <ptb...@gm...> - 2013-09-21 03:26:55
|
Hi Manav, On Fri, Sep 20, 2013 at 10:22 PM, Manav Bhatia <bha...@gm...>wrote: > Did you try this in 2D or 3D? > 2D > I do not get this error in 2D, but it shows up in 3D. > Hmm, OK. > > The following code is able to reproduce the issue at my end. > > int main (int argc, char* const argv[]) > { > // Initialize libMesh. > LibMeshInit init (argc, argv); > > > { > ParallelMesh m(init.comm()); > MeshTools::Generation::build_cube(m, 5,5,5, 0., 1., 0., 1., 0., > 1., HEX8); > EquationSystems es(m); > ExplicitSystem& s = es.add_system<ExplicitSystem>("mys"); > FEType fe; > s.add_variable("v", fe); > es.init(); > class Func: public FEMFunctionBase<Real> { > public: > virtual AutoPtr<FEMFunctionBase<Real> > clone () const { > return AutoPtr<FEMFunctionBase<Real> > (new Func); }; > virtual Real operator() (const FEMContext&, const Point& p, > const Real time = 0.) { return 0.; }; > virtual void operator() (const FEMContext&, const Point& p, > const Real time, > DenseVector<Real>& output) { > output.zero(); }; > > > }; > Func ff; > s.project_solution(&ff); > } > return 0; > } > > > Please let me know if you see the same behavior. > Thanks for the test case! I'll play with this in a little while and see if I can reproduce. Best, Paul |
From: Paul T. B. <ptb...@gm...> - 2013-09-21 06:23:37
|
On Fri, Sep 20, 2013 at 10:26 PM, Paul T. Bauman <ptb...@gm...> wrote: > > Thanks for the test case! I'll play with this in a little while and see if > I can reproduce. > OK, I can reproduce. Will figure out the problem and get a PR in. You said this worked previously - any chance you have a hash or version number associated with that? Thanks for reporting. Best, Paul |
From: Manav B. <bha...@gm...> - 2013-09-21 12:49:44
|
Hi Paul, Now that I think of it, I had never before tested this for 3d cases. Sorry that I made the wrong claim. Manav > On Sep 21, 2013, at 2:23 AM, "Paul T. Bauman" <ptb...@gm...> wrote: > > On Fri, Sep 20, 2013 at 10:26 PM, Paul T. Bauman <ptb...@gm...> wrote: >> >> >> Thanks for the test case! I'll play with this in a little while and see if I can reproduce. > > OK, I can reproduce. Will figure out the problem and get a PR in. You said this worked previously - any chance you have a hash or version number associated with that? > > Thanks for reporting. > > Best, > > Paul |
From: Paul T. B. <ptb...@gm...> - 2013-09-21 15:17:12
|
No problem! This is still good info - probably means this has been broken since I originally wrote that code. On Sep 21, 2013, at 7:49 AM, Manav Bhatia <bha...@gm...> wrote: > Hi Paul, > > Now that I think of it, I had never before tested this for 3d cases. > > Sorry that I made the wrong claim. > > Manav > > > > On Sep 21, 2013, at 2:23 AM, "Paul T. Bauman" <ptb...@gm...> wrote: > >> On Fri, Sep 20, 2013 at 10:26 PM, Paul T. Bauman <ptb...@gm...> wrote: >> >> Thanks for the test case! I'll play with this in a little while and see if I can reproduce. >> >> OK, I can reproduce. Will figure out the problem and get a PR in. You said this worked previously - any chance you have a hash or version number associated with that? >> >> Thanks for reporting. >> >> Best, >> >> Paul |
From: Paul T. B. <ptb...@gm...> - 2013-09-21 17:50:32
|
OK, I've got it. PR forthcoming. Manav was on the right track moving the side_fe_reinit around. reinit needs to be get called to properly initialize the quadrature rule. However, the reason this didn't fail in 2D or on the edge in 3D was because QGauss (the default quadrature rule here) internally inits the rule if the dimension is 1D (sides in 2D, edges in 3D). Fix forthcoming. Thanks again for reporting, Manav, and sorry for the trouble. Best, Paul On Sat, Sep 21, 2013 at 10:17 AM, Paul T. Bauman <ptb...@gm...> wrote: > No problem! This is still good info - probably means this has been broken > since I originally wrote that code. > > > > On Sep 21, 2013, at 7:49 AM, Manav Bhatia <bha...@gm...> wrote: > > Hi Paul, > > Now that I think of it, I had never before tested this for 3d cases. > > Sorry that I made the wrong claim. > > Manav > > > > On Sep 21, 2013, at 2:23 AM, "Paul T. Bauman" <ptb...@gm...> wrote: > > On Fri, Sep 20, 2013 at 10:26 PM, Paul T. Bauman <ptb...@gm...>wrote: > >> >> Thanks for the test case! I'll play with this in a little while and see >> if I can reproduce. >> > > OK, I can reproduce. Will figure out the problem and get a PR in. You said > this worked previously - any chance you have a hash or version number > associated with that? > > Thanks for reporting. > > Best, > > Paul > > |