libMesh: A C++ Finite Element Library

Tim,

I am currently using MeshFunction... and I did run into the problems
you spoke about... and I've worked around each in similar ways.  It's
getting pretty late now... so I will respond tomorrow with more
details, but I just wanted to get a response out there.

To everyone else: I was going to report these problems at some point,
but at the time I was knee deep in my own problems and was just trying
to get stuff working... hopefully after my comments tomorrow we can
get these small problems resolved.

Derek

On 7/28/06, Tim Kr=F6ger <tim@...> wrote:
> Dear all,
>
> During my work with libMesh, I encountered the following three
> problems with the MeshFunction class:
>
> 1.) The MeshFunction class cannot be instantiated since the clear()
> function is abstract.  My workaround: I use a derived class
> MyMeshFunction in which I implement a clear() function that does
> nothing.
>
> 2.) The constructor that takes "unsigend int var" (rather than
> "std::vector<unsigend int> vars") initializes the member _system_vars
> with "std::vector<unsigned int>(var)" which is a vector of var
> components, all being zero, rather than "std::vector<unsigned
> int>(1,var)".  My workaround is to use the other constructor.
>
> 3.) When calling the init() function, I get the message "ERROR:
> Already initialized!  Will ignore this call...", although I definitely
> call this function only once.
>
> Can anyone comment on this?
>
> Best Regards,
>
> Tim
>
> -------------------------------------------------------------------------
> Take Surveys. Earn Cash. Influence the Future of IT
> Join SourceForge.net's Techsay panel and you'll get the chance to share y=
our
> opinions on IT & business topics through brief surveys -- and earn cash
> http://www.techsay.com/default.php?page=3Djoin.php&p=3Dsourceforge&CID=3D=
DEVDEV
> _______________________________________________
> Libmesh-users mailing list
> Libmesh-users@...
> https://lists.sourceforge.net/lists/listinfo/libmesh-users
>

> 1.) The MeshFunction class cannot be instantiated since the clear()
> function is abstract.  My workaround: I use a derived class
> MyMeshFunction in which I implement a clear() function that does
> nothing.

This is essentially what I did... except I just added it to the actual
class itself... in the header file.  I've attached the diff.

>
> 2.) The constructor that takes "unsigend int var" (rather than
> "std::vector<unsigend int> vars") initializes the member _system_vars
> with "std::vector<unsigned int>(var)" which is a vector of var
> components, all being zero, rather than "std::vector<unsigned
> int>(1,var)".  My workaround is to use the other constructor.

I don't know about this... I must have used the other constructor all
along.  This is probably something to be looked into.

>
> 3.) When calling the init() function, I get the message "ERROR:
> Already initialized!  Will ignore this call...", although I definitely
> call this function only once.

All you have to do is comment out this->_point_locator->init(); in
mesh_function.C.  It looks like the point_locator class does this
automatically now... so I think this is crusty code that wasn't
updated when point_locator was... the problem is it is trying to
initialize it twice...  I've attached my diff for that as well.

> Best Regards,
>
> Tim

Hope that helps!

Also... to everyone else.... why is it that mesh_funciton.C is in
numerics when mesh_function.h is in mesh?  I personally vote for both
being in mesh!

Derek

Thanks, Derek, for your reply.

Could please someone who has write access to the CVS repository check 
your changes in, so that other peeople don't run into the same 
problems?

On Sat, 29 Jul 2006, Derek Gaston wrote:

>> 2.) The constructor that takes "unsigend int var" (rather than
>> "std::vector<unsigend int> vars") initializes the member _system_vars
>> with "std::vector<unsigned int>(var)" which is a vector of var
>> components, all being zero, rather than "std::vector<unsigned
>> int>(1,var)".  My workaround is to use the other constructor.
>
> I don't know about this... I must have used the other constructor all
> along.  This is probably something to be looked into.

It's a very easy task: In line 65 of mesh_function.C, the "(var)" has 
to be changed into "(1,var)".  Could someone please also check in this 
change?

> Also... to everyone else.... why is it that mesh_funciton.C is in
> numerics when mesh_function.h is in mesh?  I personally vote for both
> being in mesh!

I agree.

Best Regards,

Tim

Dear all,

Does the MeshFunction class work correctly in parallel?  I.e., do all 
processors get the correct result for the value/gradient at a given 
point or only the processor on which the element lives that contains 
the point?

Best Regards,

Tim

-- 
Dr. Tim Kroeger                            Phone +49-421-218-7710
tim@..., tim@...      Fax   +49-421-218-4236

MeVis Research GmbH, Universitaetsallee 29, D-28359 Bremen, Germany

Amtsgericht Bremen HRB 16222
Geschaeftsfuehrer: Prof. Dr. H.-O. Peitgen

On Thu, 5 Jul 2007, Tim Kröger wrote:

> Does the MeshFunction class work correctly in parallel?  I.e., do all
> processors get the correct result for the value/gradient at a given
> point or only the processor on which the element lives that contains
> the point?

Currently that's entirely up to you.  If you give MeshFunction a
standard distributed vector then you'd better only call it with points
that live on the local processor.  If you need to evaluate points that
live on any processor, you need to give MeshFunction a serial vector
to use, doing any necessary NumericVector::localize() work manually.
---
Roy

Dear all,

At the end of my application, I want to evaluate the computed solution 
at a number of points (which cannot be assumed to be grid nodes). 
Since I only want to write the values to a file, it suffices to do 
this on one processor.

My idea was:

1. Localize.
2. On all processors except #0, let everything go out of scope.
3. On processor #0, make a MeshFunction and evaluate as required.

However, this fails to work because there is a call to parallel_only() 
in the initialization of MeshFunction (that is, in tree.C, line 45, 
MeshTools::bounding_box() is called, which itself calls 
parallel_only()).

Is there any chance to work around this problem?

An idea would be to allow the user to manually specify the bounding 
box of the mesh.

Any suggestions welcome.

Best Regards,

Tim

-- 
Dr. Tim Kroeger
tim.kroeger@...            Phone +49-421-218-7710
tim.kroeger@...            Fax   +49-421-218-4236

Fraunhofer MEVIS, Institute for Medical Image Computing
Universitaetsallee 29, 28359 Bremen, Germany

On Mon, 12 Jan 2009, Tim Kroeger wrote:

> At the end of my application, I want to evaluate the computed solution
> at a number of points (which cannot be assumed to be grid nodes).
> Since I only want to write the values to a file, it suffices to do
> this on one processor.
>
> My idea was:
>
> 1. Localize.
> 2. On all processors except #0, let everything go out of scope.

Any reason why?  If you're not needing the memory for anything else,
it seems like "don't let everything go out of scope, just create a
MeshFunction or call MeshTools::bounding_box()" would be the simplest
workaround.

> 3. On processor #0, make a MeshFunction and evaluate as required.

If you build a MeshFunction but only intend to use it on proc 0,
there's no need for the only expensive step of serializing the
solution vector to the rest of those processors; I doubt you'd even
have to allocate a zeroed serial vector on them.

> However, this fails to work because there is a call to
> parallel_only() in the initialization of MeshFunction (that is, in
> tree.C, line 45, MeshTools::bounding_box() is called, which itself
> calls parallel_only()).

Right - in the case of a ParallelMesh, for instance, it's not even
possible to compute a bounding box without every processor's
assistance.

Really, what we need is to properly parallelize MeshFunction, but I'm
not sure what the right API is to do that with MPI.  If processor A
wants to evaluate the function on a point in processor B's elements,
processor B (which probably has a whole list of evaluations it needs
itself) is going to need to be notified somehow to interrupt what it's
doing and help A.
---
Roy

On Jan 12, 2009, at 10:05 AM, Roy Stogner wrote:

> Really, what we need is to properly parallelize MeshFunction, but I'm
> not sure what the right API is to do that with MPI.  If processor A
> wants to evaluate the function on a point in processor B's elements,
> processor B (which probably has a whole list of evaluations it needs
> itself) is going to need to be notified somehow to interrupt what it's
> doing and help A.

We had this capability in Sierra... the algorithm went something like  
this:

1.  Build bounding boxes on all processors that encompasses the  
elements held by that processor.
2.  Each processor communicates to every other processor it's bounding  
box.
3.  Each processor does a coarse search using the bounding boxes to  
group the points it's interested in by processor.
4.  The points most likely to fall on each processor are transmitted  
to each processor.
5.  A fine grained search on every processor is done to insure that  
the points actually do fall in an element that processor owns.  If it  
doesn't then a second guess is made at what processor it fits on and  
another round of communication happens.
6.  After every processor has the set of points to be evaluated that  
fall in it's elements... it then does all of the evaluations.
7.  The evaluations are then broadcast back to the original processor  
that needed the information.

Yes.... this is fairly involved.... and there might be a better way.   
There was quite a bit of code required to make all of this work (as  
you might imagine... this is Sierra of course!)... but it _did_ work.   
We used this capability in Encore frequently for our postprocessing  
needs.

Derek

On Mon, 12 Jan 2009, Derek Gaston wrote:

> On Jan 12, 2009, at 10:05 AM, Roy Stogner wrote:
>
>> Really, what we need is to properly parallelize MeshFunction, but I'm
>> not sure what the right API is to do that with MPI.  If processor A
>> wants to evaluate the function on a point in processor B's elements,
>> processor B (which probably has a whole list of evaluations it needs
>> itself) is going to need to be notified somehow to interrupt what it's
>> doing and help A.
>
> We had this capability in Sierra... the algorithm went something like this:
>
> 1.  Build bounding boxes on all processors that encompasses the elements held 
> by that processor.
> 2.  Each processor communicates to every other processor it's bounding box.
> 3.  Each processor does a coarse search using the bounding boxes to group the 
> points it's interested in by processor.
> 4.  The points most likely to fall on each processor are transmitted to each 
> processor.
> 5.  A fine grained search on every processor is done to insure that the 
> points actually do fall in an element that processor owns.  If it doesn't 
> then a second guess is made at what processor it fits on and another round of 
> communication happens.
> 6.  After every processor has the set of points to be evaluated that fall in 
> it's elements... it then does all of the evaluations.
> 7.  The evaluations are then broadcast back to the original processor that 
> needed the information.
>
> Yes.... this is fairly involved.... and there might be a better way.

Actually, I suspect this is the best way to do things (with MPI, at
least).  And the amount of work to be done for it isn't too scary
after writing the ParallelMesh synchronization code.  What's held me
off from writing it is that it's not a parallelization of the old API,
it would require a wholly new API, so it's not a drop-in fix for old
code.

On the other hand, maybe that's a reason to write it now: put in the
new API, write the quick serial implementation for it, and mark the
individual evaluation function as deprecated() to encourage users to
switch.  Then write the parallel implementation later.
---
Roy

Dear Roy,

On Mon, 12 Jan 2009, Roy Stogner wrote:

> On Mon, 12 Jan 2009, Tim Kroeger wrote:
>
>> At the end of my application, I want to evaluate the computed solution
>> at a number of points (which cannot be assumed to be grid nodes).
>> Since I only want to write the values to a file, it suffices to do
>> this on one processor.
>> 
>> My idea was:
>> 
>> 1. Localize.
>> 2. On all processors except #0, let everything go out of scope.
>
> Any reason why?  If you're not needing the memory for anything else,

Well, I thought that creating the MeshFunction would require a lot of 
memory, so it would be nice to let the other processors reduce their 
memory requirements.  (Some of the processors are located on the same 
node, i.e. are using the same memory pool.)

> it seems like "don't let everything go out of scope, just create a
> MeshFunction or call MeshTools::bounding_box()" would be the simplest
> workaround.

That's what I'm doing now, and it works, but I thought there would be 
a nicer solution.

>> 3. On processor #0, make a MeshFunction and evaluate as required.
>
> If you build a MeshFunction but only intend to use it on proc 0,
> there's no need for the only expensive step of serializing the
> solution vector to the rest of those processors; I doubt you'd even
> have to allocate a zeroed serial vector on them.

Well, the point is that there is no NumericVector::localize_to_one 
(NumericVector&, unsigned int).  See also my posting of December 17, 
2008.

>> However, this fails to work because there is a call to
>> parallel_only() in the initialization of MeshFunction (that is, in
>> tree.C, line 45, MeshTools::bounding_box() is called, which itself
>> calls parallel_only()).
>
> Right - in the case of a ParallelMesh, for instance, it's not even
> possible to compute a bounding box without every processor's
> assistance.

That's true, of course.

Best Regards,

Tim

-- 
Dr. Tim Kroeger
tim.kroeger@...            Phone +49-421-218-7710
tim.kroeger@...            Fax   +49-421-218-4236

Fraunhofer MEVIS, Institute for Medical Image Computing
Universitaetsallee 29, 28359 Bremen, Germany

Dear all,

Can someone please briefly explain how to use MeshFunction to me? I've been digging the mailing list archives and ./src/error_estimation/exact_error_estimator.C to understand it but I'm still hazy on how it use it? 

This is what I intend to do I want to get the L-2 projection of a solution on a side of a 3-D mesh on 2-D mesh. 

PS: feel free to comment if you think there is better way to go about doing the L-2 projection from 3-D to 2-D.

Best,
Ata

Hi, 

  Does MeshFunction expect the vector (input argument for constructor) to be serialized and hold the global solution? 

  I am attempting to use the class in my code. Passing it the system.solution without serializing it is leading to nans as the interpolated solution at the boundary of each partition on processor_id != 0. 

  I looked at the code in exact_solution.C and exact_error_estimator.C, both of which serialize the solution before passing it to MeshFunction. Not sure why this would be necessary. 

Thanks, 
Manav

Yes.

Ok, well, MeshFunction itself doesn't care (as you found out)... but
if you want to evaluate your solution in elements containing
off-processor degrees of freedom then you need to serialize that
vector.

Think of it this way: MeshFunction doesn't do parallel
communication... you are responsible for feeding it the values you
need.

Derek

Sent from my iPad

On May 24, 2013, at 7:54 AM, Manav Bhatia <bhatiamanav@...> wrote:

> Hi,
>
>  Does MeshFunction expect the vector (input argument for constructor) to be serialized and hold the global solution?
>
>  I am attempting to use the class in my code. Passing it the system.solution without serializing it is leading to nans as the interpolated solution at the boundary of each partition on processor_id != 0.
>
>  I looked at the code in exact_solution.C and exact_error_estimator.C, both of which serialize the solution before passing it to MeshFunction. Not sure why this would be necessary.
>
> Thanks,
> Manav
>
>
>
> ------------------------------------------------------------------------------
> Try New Relic Now & We'll Send You this Cool Shirt
> New Relic is the only SaaS-based application performance monitoring service
> that delivers powerful full stack analytics. Optimize and monitor your
> browser, app, & servers with just a few lines of code. Try New Relic
> and get this awesome Nerd Life shirt! http://p.sf.net/sfu/newrelic_d2d_may
> _______________________________________________
> Libmesh-users mailing list
> Libmesh-users@...
> https://lists.sourceforge.net/lists/listinfo/libmesh-users

On May 24, 2013, at 9:58 AM, Derek Gaston <friedmud@...> wrote:

> Yes.
> 
> Ok, well, MeshFunction itself doesn't care (as you found out)... but
> if you want to evaluate your solution in elements containing
> off-processor degrees of freedom then you need to serialize that
> vector.
> 


But if the element is on a given processor, then wouldn't the local portion of system.solution contain all the relevant dofs to calculate the solution on that element? 

If so, then passing it the system.solution should be alright. No? 


> Think of it this way: MeshFunction doesn't do parallel
> communication... you are responsible for feeding it the values you
> need.
> 
> Derek
> 
> Sent from my iPad
>

What are you trying to do?  There are many scenarios where this could
fail... and some where it might succeed.

Derek

Sent from my iPad

On May 24, 2013, at 8:02 AM, Manav Bhatia <bhatiamanav@...> wrote:

>
>
> On May 24, 2013, at 9:58 AM, Derek Gaston <friedmud@...> wrote:
>
>> Yes.
>>
>> Ok, well, MeshFunction itself doesn't care (as you found out)... but
>> if you want to evaluate your solution in elements containing
>> off-processor degrees of freedom then you need to serialize that
>> vector.
>
>
> But if the element is on a given processor, then wouldn't the local portion of system.solution contain all the relevant dofs to calculate the solution on that element?
>
> If so, then passing it the system.solution should be alright. No?
>
>
>> Think of it this way: MeshFunction doesn't do parallel
>> communication... you are responsible for feeding it the values you
>> need.
>>
>> Derek
>>
>> Sent from my iPad
>

Put another way, solution contains storage for remote entries but it cannot keep them in sync under all usage cases (1) for performance reasons and (2) to avoid parallel deadlocks. 

So depending on when you are accessing it the remote values may be stale, leading to chaos. 

-Ben



On May 24, 2013, at 9:08 AM, "Derek Gaston" <friedmud@...> wrote:

> What are you trying to do?  There are many scenarios where this could
> fail... and some where it might succeed.
> 
> Derek
> 
> Sent from my iPad
> 
> On May 24, 2013, at 8:02 AM, Manav Bhatia <bhatiamanav@...> wrote:
> 
>> 
>> 
>> On May 24, 2013, at 9:58 AM, Derek Gaston <friedmud@...> wrote:
>> 
>>> Yes.
>>> 
>>> Ok, well, MeshFunction itself doesn't care (as you found out)... but
>>> if you want to evaluate your solution in elements containing
>>> off-processor degrees of freedom then you need to serialize that
>>> vector.
>> 
>> 
>> But if the element is on a given processor, then wouldn't the local portion of system.solution contain all the relevant dofs to calculate the solution on that element?
>> 
>> If so, then passing it the system.solution should be alright. No?
>> 
>> 
>>> Think of it this way: MeshFunction doesn't do parallel
>>> communication... you are responsible for feeding it the values you
>>> need.
>>> 
>>> Derek
>>> 
>>> Sent from my iPad
> 
> ------------------------------------------------------------------------------
> Try New Relic Now & We'll Send You this Cool Shirt
> New Relic is the only SaaS-based application performance monitoring service 
> that delivers powerful full stack analytics. Optimize and monitor your
> browser, app, & servers with just a few lines of code. Try New Relic
> and get this awesome Nerd Life shirt! http://p.sf.net/sfu/newrelic_d2d_may
> _______________________________________________
> Libmesh-users mailing list
> Libmesh-users@...
> https://lists.sourceforge.net/lists/listinfo/libmesh-users

I tried running in the debug mode (should have done that earlier) and the code aborted with the following error. 


No index 128 in ghosted vector.
Vector contains [1352,2604)
And empty ghost array.
[1] ./include/libmesh/petsc_vector.h, line 1176, compiled May 22 2013 at 22:27:31


It also was trying to get value for an off-processor element. So, this perhaps explains the hans for the optimized executable. 

I was assuming that the project solution will work with only the local elements, but that does not seem to be the case. Will try with serialization of the vector. 

Manav




On May 24, 2013, at 10:08 AM, Derek Gaston <friedmud@...> wrote:

> What are you trying to do?  There are many scenarios where this could
> fail... and some where it might succeed.
> 
> Derek
> 
> Sent from my iPad
> 
> On May 24, 2013, at 8:02 AM, Manav Bhatia <bhatiamanav@...> wrote:
> 
>> 
>> 
>> On May 24, 2013, at 9:58 AM, Derek Gaston <friedmud@...> wrote:
>> 
>>> Yes.
>>> 
>>> Ok, well, MeshFunction itself doesn't care (as you found out)... but
>>> if you want to evaluate your solution in elements containing
>>> off-processor degrees of freedom then you need to serialize that
>>> vector.
>> 
>> 
>> But if the element is on a given processor, then wouldn't the local portion of system.solution contain all the relevant dofs to calculate the solution on that element?
>> 
>> If so, then passing it the system.solution should be alright. No?
>> 
>> 
>>> Think of it this way: MeshFunction doesn't do parallel
>>> communication... you are responsible for feeding it the values you
>>> need.
>>> 
>>> Derek
>>> 
>>> Sent from my iPad
>>