libmesh-devel — General information for libmesh developers

[Libmesh-devel] Profiling libmesh

[Fabio C. <fab...@gm...>]

Hey experts =),

I`m going some profiling of our code with focus on energy and just wanted
to know if there is any recommended method, library or something that has
been a wish in this regards.

Regards,
Fabio C. Canesin

Re: [Libmesh-devel] Multigrid Projection Matrixes with Periodic Boundaries

[Roy S. <roy...@ic...>]

On Thu, 6 Jul 2017, Boris Boutkov wrote:

> Hello again - Ive been doing further work on the multigrid
> projection matrix front and things are starting to look fairly
> stable overall!

Great!

> The main issue that has resurfaced is the treatment of the boundary
> conditions. Though it seems that for most cases convergence rates
> are as expected, there clearly still remains a problem with the
> periodic cases in even the simplest 1D/2D test problems.

Strange.  I suppose I can imagine 2D giving us trouble, but in 1D,
every restriction or prolongation of a periodic solution will still be
a periodic solution.

On the other hand...  In the basic multigrid methods, we do the
prolongation on a coarse error term, and the restriction on a fine
residual.

The error term should be periodic, so we shouldn't need to care if the
prolongation doesn't redundantly enforce that.

The fine residual, though... it's actually hard for me to picture in
my head what prolongation really means in the PDE case, because
although the residual coefficients have the same algebraic structure
as the primal solution coefficients, the residual is in the dual
space.  I don't think "pretend they're coefficients to a function and
project that" is the correct thing to do, but I worry that maybe it's
correct enough to work in simple cases and incorrect enough to fail as
we add complications like periodicity and AMR.  "Find the function
which gives the same L2 inner product as the residual vector's l2
inner product, then project that function, then find the residual
vector in the new space which gives the same l2 inner product as the
projected function's L2 inner product", like a discrete Riesz
representation transformation, maybe?

Anyway, that's an aside.  However the residual vector is being
restricted, I take it it's not already conforming to the periodic
constraints, even in 1D, and so we need to modify the matrix to make
sure it's conforming after the restriction?

> I think that the natural way forward is that we need to constrain
> the element matrixes as we are building up the projection matrixes,

Well, we may need an option to do so, at least.  It definitely needs
to be a functor that we can make an inline no-op in the common cases,
lest that refactored projection be slowed down again.

> the first place that comes to mind would be right after we have set
> the target_matrix dof information inside the MatrixFillAction
> component of the GenericProjector.

I'm not sure if this is right, but it's at least probably worth trying.

> The issue then is that the usual DofMap::constrain_element_matrix
> implementations only have DenseMatrix versions whereas here it seems
> we would need Sparse ones in the projection matrix situation.

True.  Can we get away with the same "move the implementations to a
templated helper function, then instantiate the common cases in the
original .C file for backwards compatibility, then instantiate the
DSNA cases in the new projection matrix code" trick?
---
Roy

[Libmesh-devel] Multigrid Projection Matrixes with Periodic Boundaries

[Boris B. <bor...@bu...>]

Hello again - Ive been doing further work on the multigrid projection
matrix front and things are starting to look fairly stable overall!

The main issue that has resurfaced is the treatment of the boundary
conditions. Though it seems that for most cases convergence rates are as
expected, there clearly still remains a problem with the periodic cases in
even the simplest 1D/2D test problems. I think that the natural way forward
is that we need to constrain the element matrixes as we are building up the
projection matrixes, the first place that comes to mind would be right
after we have set the target_matrix dof information inside the
MatrixFillAction component of the GenericProjector. The issue then is that
the usual  DofMap::constrain_element_matrix implementations only have
DenseMatrix versions whereas here it seems we would need Sparse ones in the
projection matrix situation.

In essence I'd like to verify that creating these analagous DofMap calls is
the right track to follow as Im still a bit shakey with the BC and DSNA
components of the library - so I'd appreciate any input as to the viability
of this gameplan or if an altogether different course of action is needed.


Thanks in advance for the input,
Boris

Re: [Libmesh-devel] Raviart Thomas

[Roy S. <roy...@ic...>]

On Mon, 5 Jun 2017, Rovi Gabriele wrote:

> I would like to implement the Raviart Thomas elements. In order to obtain a global basis function, if the global
> dofs of the actual element are such that:
> 
> n1 < n2 < n3 (2D) or n1 < n2 < n3 < n4 (3D)
> 
> then we have the mapping between:
> 
>  0 <——> n1
>  1 <——> n2
>  2 <——> n3
> (3 <——> n4, 3D)
> 
> This kind of mapping can have a positive or negative determinant. So it can be that a clockwise element is mapped in
> the reference counterclockwise one.

In libMesh we require 3D elements (as well as 2D elements, if the
library is configured --2d-only and so we can rule out 2D manifolds
curving through R3) to have positive determinants.

> Where is defined (header and source with the computation) the mapping between the reference and actual elements? 

Elem::node_id(i) (in elem.h) takes local node id i (0,1,2,3 above) and
returns the global node id (n1, n2, n3, n4 above).  There is never a
guarantee about the ordering of node ids, however.  To obtain globally
unique orderings of edges and faces for FE types which require that,
the current procedure is: compare the nodes physical locations, and
hope the user isn't using mesh motion on this problem.  Search for
"edge-flipping" in fe_hierarchic_shape_2D.C for the simplest example.

Obviously "hope the user isn't using mesh motion on this problem"
is a bug to be fixed here; the "unique_id()" facility in libMesh is
probably our best bet to change that, but that's an optional feature
at the moment, and because it noticeably adds to memory usage it will
probably remain optional-only indefinitely.

> Is it defined as “counterclockwise <—>counterclockwise? If so, would be a drastic change to define another one like
> the one that I have now described?

Trying to force any particular node numbering would be a very drastic
change; we do have a "don't renumber" option for users who expect
their original mesh numbering to remain unchanged through adaptivity,
but even in that case newly created nodes have to be given new
available numberings (i.e. probably higher numbers than have been
already assigned) which may conflict with your desired ordering.
---
Roy

Re: [Libmesh-devel] OpenHPC

[Julian A. <ju...@tf...>]

On Fri, Jun 23, 2017 at 5:55 PM, Roy Stogner <roy...@ic...> wrote:
>
> On Thu, 22 Jun 2017, Julian Andrej wrote:
>
>> Isn't libmesh a component which is heavily dependent on the users
>> needs? So a "one size fits all" RPM is a bit off.
>
>
> For most of our configuration options, the relevant users' needs are
> "I need to not spend a lot of time hunting down and reading about and
> configuring and building a dozen other packages first".  This need is
> already satisfied, even with a fully-featured libMesh, if "yum install
> libMesh" automatically installs all those features' dependencies too.
>
> We'll want "two sizes fit all", IMHO - "libmesh" and
> "libmesh-complex".

I agree with that after thinking once more about the "general user".

>
> (that said, Spack looks cool)
> ---
> Roy



-- 
Julian Andrej

Chair of Automatic Control
Faculty of Engineering
Kiel University
Kaiserstrasse 2 | 24143 Kiel | Germany

L: Room F-116
T: +49(0)431 880-6121
F: +49(0)431 880-6278
ju...@tf... | http://www.control.tf.uni-kiel.de

Re: [Libmesh-devel] OpenHPC

[Roy S. <roy...@ic...>]

On Thu, 22 Jun 2017, Julian Andrej wrote:

> Isn't libmesh a component which is heavily dependent on the users
> needs? So a "one size fits all" RPM is a bit off.

For most of our configuration options, the relevant users' needs are
"I need to not spend a lot of time hunting down and reading about and
configuring and building a dozen other packages first".  This need is
already satisfied, even with a fully-featured libMesh, if "yum install
libMesh" automatically installs all those features' dependencies too.

We'll want "two sizes fit all", IMHO - "libmesh" and
"libmesh-complex".

(that said, Spack looks cool)
---
Roy

Re: [Libmesh-devel] OpenHPC

[Kirk, B. (JSC-EG311) <ben...@na...>]

Thanks for the link!

I really don’t see it as an either/or; both would certainly be welcome and are valuable for their own audiences.

As for the premise of the question, I’m not sure.  Certainly libMesh supports a lot of configure options (too many now?), but a default build compiled against a binary package of PETSc would be pretty easy to do and valuable for someone just getting in to the library.

-Ben


> On Jun 22, 2017, at 2:30 PM, Julian Andrej <ju...@tf...> wrote:
> 
> Isn't libmesh a component which is heavily dependent on the users
> needs? So a "one size fits all" RPM is a bit off.
> 
> In my opinion extending something like [1] (to support more options)
> would be more appropriate. With Spack you can choose whatever
> compiler/MPI/PETSc/Trilinos version or version combination YOU want to
> have.
> 
> [1] https://github.com/LLNL/spack/blob/develop/var/spack/repos/builtin/packages/libmesh/package.py
> 
> On Thu, Jun 22, 2017 at 9:21 PM, John Peterson <jwp...@gm...> wrote:

Re: [Libmesh-devel] OpenHPC

[Julian A. <ju...@tf...>]

Isn't libmesh a component which is heavily dependent on the users
needs? So a "one size fits all" RPM is a bit off.

In my opinion extending something like [1] (to support more options)
would be more appropriate. With Spack you can choose whatever
compiler/MPI/PETSc/Trilinos version or version combination YOU want to
have.

[1] https://github.com/LLNL/spack/blob/develop/var/spack/repos/builtin/packages/libmesh/package.py

On Thu, Jun 22, 2017 at 9:21 PM, John Peterson <jwp...@gm...> wrote:
>
>
> On Thu, Jun 22, 2017 at 1:02 PM, Kirk, Benjamin (JSC-EG311)
> <ben...@na...> wrote:
>>
>> Are any of you guys aware of Karl’s project:
>> http://openhpc.community/
>>
>> A nice yum-based way to kickstart and maintain an HPC environment
>> consistent with what a lot of us now like to see.
>>
>> I’m considering contributing a libmesh build using their framework, but
>> was curious if you guys have any experience running it yet.  Its on my list
>> now but might be a bit...
>
>
> That URL is... interesting. I've been wondering when I would see one of the
> new TLDs in real life.
>
> I haven't clicked around on the site much yet, but it sounds like we would
> fit well in the parallel-libs "Component".
>
> All their stuff appears to be RPM-based? I recently removed a very old
> libmesh.spec file in a0a0d81a, maybe that would be a good starting point for
> any work on this...
>
> --
> John
>
> ------------------------------------------------------------------------------
> Check out the vibrant tech community on one of the world's most
> engaging tech sites, Slashdot.org! http://sdm.link/slashdot
> _______________________________________________
> Libmesh-devel mailing list
> Lib...@li...
> https://lists.sourceforge.net/lists/listinfo/libmesh-devel
>



-- 
Julian Andrej

Chair of Automatic Control
Faculty of Engineering
Kiel University
Kaiserstrasse 2 | 24143 Kiel | Germany

L: Room F-116
T: +49(0)431 880-6121
F: +49(0)431 880-6278
ju...@tf... | http://www.control.tf.uni-kiel.de

Re: [Libmesh-devel] OpenHPC

[Kirk, B. (JSC-EG311) <ben...@na...>]

On Jun 22, 2017, at 2:21 PM, John Peterson <jwp...@gm...<mailto:jwp...@gm...>> wrote:

All their stuff appears to be RPM-based? I recently removed a very old libmesh.spec file in a0a0d81a, maybe that would be a good starting point for any work on this...


Yeah, I peeked at the petsc spec and it seems pretty simple. I think it would be an easy thing to start from.

I’ll be insalling a vanilla 7.3 machine here in the near future and may try layering this on top for MPI and supporting libraries.

Re: [Libmesh-devel] OpenHPC

[John P. <jwp...@gm...>]

On Thu, Jun 22, 2017 at 1:02 PM, Kirk, Benjamin (JSC-EG311) <
ben...@na...> wrote:

> Are any of you guys aware of Karl’s project:
> http://openhpc.community/
>
> A nice yum-based way to kickstart and maintain an HPC environment
> consistent with what a lot of us now like to see.
>
> I’m considering contributing a libmesh build using their framework, but
> was curious if you guys have any experience running it yet.  Its on my list
> now but might be a bit...
>

That URL is... interesting. I've been wondering when I would see one of the
new TLDs in real life.

I haven't clicked around on the site much yet, but it sounds like we would
fit well in the parallel-libs "Component".

All their stuff appears to be RPM-based? I recently removed a very old
libmesh.spec file in a0a0d81a, maybe that would be a good starting point
for any work on this...

-- 
John

Re: [Libmesh-devel] OpenHPC

[Derek G. <fri...@gm...>]

I haven't heard anything about it.  If you end up trying it out I would be
interested in hearing about how well it works.  There are several clusters
here at MIT where we've had to roll our own HPC distro...

Derek

On Thu, Jun 22, 2017 at 3:02 PM Kirk, Benjamin (JSC-EG311) <
ben...@na...> wrote:

> Are any of you guys aware of Karl’s project:
> http://openhpc.community/
>
> A nice yum-based way to kickstart and maintain an HPC environment
> consistent with what a lot of us now like to see.
>
> I’m considering contributing a libmesh build using their framework, but
> was curious if you guys have any experience running it yet.  Its on my list
> now but might be a bit...
>
> ------------------------------------------------------------------------------
> Check out the vibrant tech community on one of the world's most
> engaging tech sites, Slashdot.org! http://sdm.link/slashdot
> _______________________________________________
> Libmesh-devel mailing list
> Lib...@li...
> https://lists.sourceforge.net/lists/listinfo/libmesh-devel
>

[Libmesh-devel] Performance monitoring and instrumentation

[Fabio C. <fab...@gm...>]

Dear experts,

Is there anybody that used Score-P or Tau in a libMesh application ? Did
you got accurate network cost measurements ? How did you instrument,
rebuild libMesh with it ?

Thanks for any help!
Regards,
Fábio

[Libmesh-devel] Raviart Thomas

[Rovi G. <gab...@us...>]

Hi,


I would like to implement the Raviart Thomas elements. In order to obtain a global basis function, if the global dofs of the actual element are such that:


n1 < n2 < n3 (2D) or n1 < n2 < n3 < n4 (3D)


then we have the mapping between:


 0 <——> n1

 1 <——> n2

 2 <——> n3

(3 <——> n4, 3D)



This kind of mapping can have a positive or negative determinant. So it can be that a clockwise element is mapped in the reference counterclockwise one.


Where is defined (header and source with the computation) the mapping between the reference and actual elements?

Is it defined as “counterclockwise <—>counterclockwise? If so, would be a drastic change to define another one like the one that I have now described?



To be consistent, the mapping must be  unique for all the shape functions (for example, the Lagrangian ones). Does this mean that I should change also these shape functions?


How do you suggest to proceed?


Thank you in advance


GR

[Libmesh-devel] libMesh 1.2.0 released

[John P. <jwp...@gm...>]

Dear all,

LibMesh 1.2.0 has been tagged and the tarballs are available here:

https://github.com/libMesh/libmesh/releases/tag/v1.2.0

The v1.2.0 branch was initially created from master on March 17, 2017, and
several bugfixes have been cherry-picked onto it since that time. A more
detailed list of changes from 1.1.0 -> 1.2.0 is available at the link below:

https://github.com/libMesh/libmesh/blob/v1.2.0/NEWS

Thanks,
John

Re: [Libmesh-devel] [Libmesh-users] -snes_view

[David K. <dav...@ak...>]

On Wed, Apr 26, 2017 at 2:58 PM, John Peterson <jwp...@gm...> wrote:

>
>
> On Wed, Apr 26, 2017 at 12:54 PM, Fabio Canesin <fab...@gm...>
> wrote:
>
>> Dear users/developers,
>>
>> I noticed that -snes_view does not print the data it should on solvers.
>> -ksp_view works as expected.
>>
>
> Are you definitely using a SNES object (via PetscNonlinearSolver)?
>
> Not all of the libmesh examples use one...
>


FYI, I just ran systems_of_equations_ex8 (which does use
PetscNonlinearSolver via a NonlinearImplicitSystem) with -snes_view and it
showed the expected SNES output.

David

Re: [Libmesh-devel] [Libmesh-users] -snes_view

[John P. <jwp...@gm...>]

On Wed, Apr 26, 2017 at 1:05 PM, David Knezevic <dav...@ak...>
wrote:

> On Wed, Apr 26, 2017 at 2:58 PM, John Peterson <jwp...@gm...>
> wrote:
>
>>
>>
>> On Wed, Apr 26, 2017 at 12:54 PM, Fabio Canesin <fab...@gm...>
>> wrote:
>>
>>> Dear users/developers,
>>>
>>> I noticed that -snes_view does not print the data it should on solvers.
>>> -ksp_view works as expected.
>>>
>>
>> Are you definitely using a SNES object (via PetscNonlinearSolver)?
>>
>> Not all of the libmesh examples use one...
>>
>
>
> FYI, I just ran systems_of_equations_ex8 (which does use
> PetscNonlinearSolver via a NonlinearImplicitSystem) with -snes_view and it
> showed the expected SNES output.
>

Hmm... I think I agree with (both of) you, I ran adjoints_ex3 (which uses
PetscDiffSolver) with -snes_view and it didn't print anything.

I don't know of anything at the libmesh level that would prevent -snes_view
from firing though...

-- 
John

Re: [Libmesh-devel] [Libmesh-users] -snes_view

[John P. <jwp...@gm...>]

On Wed, Apr 26, 2017 at 12:54 PM, Fabio Canesin <fab...@gm...>
wrote:

> Dear users/developers,
>
> I noticed that -snes_view does not print the data it should on solvers.
> -ksp_view works as expected.
>

Are you definitely using a SNES object (via PetscNonlinearSolver)?

Not all of the libmesh examples use one...

-- 
John

[Libmesh-devel] -snes_view

[Fabio C. <fab...@gm...>]

Dear users/developers,

I noticed that -snes_view does not print the data it should on solvers.
-ksp_view works as expected.

Is this expected ?

Regards,
Fábio

Re: [Libmesh-devel] Implementing multigrid style projections through GenericProjector interface

[Boris B. <bor...@bu...>]

Awesome thanks for the updates. Ill pull in the changes and give them a try
in the next couple of days - will post back here if there are any major
hicups.

Thanks for the changes & clarifications,
- Boris

On Fri, Mar 24, 2017 at 6:23 PM, Roy Stogner <roy...@ic...>
wrote:

>
> On Wed, 22 Mar 2017, Roy Stogner wrote:
>
> On Mon, 20 Mar 2017, Roy Stogner wrote:
>>
>>  I'm starting a branch to do some of these things; I'll let you know
>>>  after I push to GitHub.
>>>
>>
>> You can see the start of things at roystgnr/projection_matrix
>>
>> Some questions that are coming up:
>>
>> In the complex-valued problem case, are we going to want projection
>> matrices to be SparseMatrix<Number> or can we still get away with
>> SparseMatrix<Real>?
>>
>> In the heterogeneous Dirichlet boundary condition case, prolongation
>> operations are *not* linear operators unless the Dirichlet BC is
>> contained in the trace of the finite element space.  I'm going to
>> cross my fingers and hope we can ignore that in the context of
>> multigrid solvers, but if anyone knows of any papers that discuss this
>> issue I'd love to hear about it.
>>
>
> And some more questions, this time directed at my senior libMesh
> developers:
>
> Whose idea was SparseMatrix::attach_dof_map??  Were we *trying* to
> stop people from ever creating more than just the one sort of square
> matrix?
>
> I finished the code for filling our projection matrices, I dug into
> the APIs to figure out what I needed to do to get the sparsity pattern
> right first... and I find out it's impossible with the existing API.
>
> Boris:
>
> In the next few months, I'll try to find time at some point to clean
> up the SparseMatrix/PetscMatrix interface.  You won't want to wait for
> that.
>
> In the near term, see if roystgnr/projection_matrix is enough to get
> you started, after you replace each of the APIs-which-use-SparseMatrix
> with APIs-which-use-raw-PETSc-Mat.  I haven't tested the code for
> obvious reasons, but it's all there and all the metaprogramming
> hurdles have been jumped.
> ---
> Roy
>

Re: [Libmesh-devel] Implementing multigrid style projections through GenericProjector interface

[Roy S. <roy...@ic...>]

On Wed, 22 Mar 2017, Roy Stogner wrote:

> On Mon, 20 Mar 2017, Roy Stogner wrote:
>
>>  I'm starting a branch to do some of these things; I'll let you know
>>  after I push to GitHub.
>
> You can see the start of things at roystgnr/projection_matrix
>
> Some questions that are coming up:
>
> In the complex-valued problem case, are we going to want projection
> matrices to be SparseMatrix<Number> or can we still get away with
> SparseMatrix<Real>?
>
> In the heterogeneous Dirichlet boundary condition case, prolongation
> operations are *not* linear operators unless the Dirichlet BC is
> contained in the trace of the finite element space.  I'm going to
> cross my fingers and hope we can ignore that in the context of
> multigrid solvers, but if anyone knows of any papers that discuss this
> issue I'd love to hear about it.

And some more questions, this time directed at my senior libMesh
developers:

Whose idea was SparseMatrix::attach_dof_map??  Were we *trying* to
stop people from ever creating more than just the one sort of square
matrix?

I finished the code for filling our projection matrices, I dug into
the APIs to figure out what I needed to do to get the sparsity pattern
right first... and I find out it's impossible with the existing API.

Boris:

In the next few months, I'll try to find time at some point to clean
up the SparseMatrix/PetscMatrix interface.  You won't want to wait for
that.

In the near term, see if roystgnr/projection_matrix is enough to get
you started, after you replace each of the APIs-which-use-SparseMatrix
with APIs-which-use-raw-PETSc-Mat.  I haven't tested the code for
obvious reasons, but it's all there and all the metaprogramming
hurdles have been jumped.
---
Roy

Re: [Libmesh-devel] Implementing multigrid style projections through GenericProjector interface

[Roy S. <roy...@ic...>]

On Mon, 20 Mar 2017, Roy Stogner wrote:

> I'm starting a branch to do some of these things; I'll let you know
> after I push to GitHub.

You can see the start of things at roystgnr/projection_matrix

Some questions that are coming up:

In the complex-valued problem case, are we going to want projection
matrices to be SparseMatrix<Number> or can we still get away with
SparseMatrix<Real>?

In the heterogeneous Dirichlet boundary condition case, prolongation
operations are *not* linear operators unless the Dirichlet BC is
contained in the trace of the finite element space.  I'm going to
cross my fingers and hope we can ignore that in the context of
multigrid solvers, but if anyone knows of any papers that discuss this
issue I'd love to hear about it.
---
Roy

Re: [Libmesh-devel] Implementing multigrid style projections through GenericProjector interface

[Roy S. <roy...@ic...>]

On Fri, 17 Mar 2017, Boris Boutkov wrote:

> I've been spending some time thinking about how to neatly implement
> multigrid projections through the existing  GenericProjector
> interface and would like some input on the following points.

Sorry for the late response; I wasn't sure I understood your
questions, so I took advantage of having Paul Bauman here in person
today to interrogate him about the problems.

> - In order to properly insert local projection matrices into the
> global projection I need some additional information in the
> ProjectionAction parameter of the GenericProjector. Specifically, I
> need access to the old dof information (as well as the existing new
> dof info) as this dictates where the local projection matrix gets
> inserted into the global one.

Right!  This is why FValue will need to be a sparse array type like
MetaPhysicL::DynamicSparseNumberArray.

> While it seemed originally that I could hijack the Fvalue parameter
> information and squeeze old dof info in there, now this appears to
> be a problem since Fvalue is also used in Fe all throughout the
> Ke/Fe solves.

In theory, that's exactly what we want!  DenseMatrix::cholesky_solve()
takes a second template type, precisely so we can do things like
using a DenseMatrix<Real> to solve a problem with DenseVector<DSNA>
data and solution.

In practice... we don't even bother to instantiate that function with
any combinations other than Real+Real, Real+Complex, and (if
--enable-complex) Complex+Complex.  So we probably need to move those
function definitions into a dense_matrix_impl.h header that we can
include elsewhere for when we want to instantiate Real+CrazyOOClass.

I'm starting a branch to do some of these things; I'll let you know
after I push to GitHub.

> - Additionally, it seems like we may need to update the
> ProjectionAction interface to take in Ke and Fe as well.

No!  As long as you can solve for Ue, making it another DSNA, then
each entry of Ue is itself a simple matrix row.  No need for Ke
afterward or Fe afterward - you get the column index for each entry
from the indices in the sparse array.
---
Roy

[Libmesh-devel] Implementing multigrid style projections through GenericProjector interface

[Boris B. <bor...@bu...>]

Hello all,

I've been spending some time thinking about how to neatly implement
multigrid projections through the existing  GenericProjector interface and
would like some input on the following points.

- In order to properly insert local projection matrices into the global
projection I need some additional information in the ProjectionAction
parameter of the GenericProjector. Specifically, I need access to the old
dof information (as well as the existing new dof info) as this dictates
where the local projection matrix gets inserted into the global one. While
it seemed originally that I could hijack the Fvalue parameter information
and squeeze old dof info in there, now this appears to be a problem since
Fvalue is also used in Fe all throughout the Ke/Fe solves.

- Additionally, it seems like we may need to update the ProjectionAction
interface to take in Ke and Fe as well. This would be needed in order to
accommodate non-nodal dofs. Ke^{-1} would for example be needed to insert
into projection matrix for non nodal dofs, and likely even for Q2 space as
well.

If this line of thinking is on the right track I would appreciate any
feedback regarding how best to deal with these issues, or any tips on how I
could avoid these problems in the first place.

Thanks in advance,
Boris Boutkov

[Libmesh-devel] EquationSystem Reinit with Refinement

[Tim A. <tra...@bu...>]

I have been looking into refactoring EquationSystem::reinit() with the
original motivation of preserving refinement flags for apps (such as GRINS)
that call refine_and_coarsen() and then es->reinit().
refine/coarsen_elements() end up getting called twice (once outside reinit,
and then again inside), and the second calls clear out the JUST_REFINED and
JUST_COARSENED flags, breaking some of my GRINS QoI code.

I created 2 simple unit tests in libMesh to check refinement flags both
when the user calls refine_and_coarsen_elements() and then es->reinit()
(refine-then-reinit), as well as when the user just sets the refinement
flags and calls es->reinit() to do the actual refinement
(reinit-and-refine). These 2 tests can be found here
<https://github.com/tradowsk/libmesh/commit/a3bdd55d24a4f7cad16efd0041caa9934fb6117f#diff-b37b7959e5e90cfb6e8e35e3cf832cd5R121>
.

However, I believe I have discovered an issue with the current
implementation of ES::reinit(). In the reinit-and-refine case (my
testReinitRefinePreserveFlags test), an assert gets tripped here
<https://github.com/libMesh/libmesh/blob/master/include/base/dof_object.h#L807>.
Additionally, I modified the GRINS AMR code to do reinit-and-refine and
tripped the same assert on the first adaptive refinement of a QoI
calculation.

I was able to find a fix for the refine-then-reinit case, which will pass
my testRefineThenReinitPreserveFlags libMesh test, as well as a few tests I
have in GRINS.

But I have been unable to find a fix for the reinit-and-refine case. I have
tried:

1) testing for dof_constraints_created==true on equation_systems.C lines
246 & 270 (passed libMesh make check, but failed a GRINS regression and AMR
test)

2) Using refine_and_coarsen_elements() instead of the current 2-step
coarsen then refine (made no difference)

3) Moving the coarsen/refine_elements() calls to the beginning of reinit()
before the new systems are added, and left distribute_dofs(),
reinit_constraints(), and get_system(i)->reinit() to the very end (fails
several libmesh tests: incorrect results on split_mesh and NULL_PTR error
in another test)

Does anyone else use the reinit-and-refine case, and if so, how did you get
it working?
I'm not very familiar with this part of the libMesh code, and I would
appreciate any feedback as to whether I am missing something or any
suggestions from more knowledgeable devs.

Thanks!

-- 
*~Tim Adowski*

Re: [Libmesh-devel] parallel mesh

[John P. <jwp...@gm...>]

This message finally appeared in the spam filters for the libmesh-devel
mailing list, but I believe it has already been discussed.

--
John



On Mon, Jan 9, 2017 at 2:57 PM, Nestola Maria Giuseppina Chiara <
mar...@us...> wrote:

> Dear all,
> we are trying  to build parallel meshes. What i need is that each
> processor has just its own mesh.
>
> We are coupling a fem code with FD code thus we need to guarantee also the
> same partition.
>
> We have already implemented a routine to do this but we have some
> questions.
>
> 1)
>
> If our setup of elements already includes elements across process
> boundaries, do we need to still call gather_neighboring_elements()?
> I added the following line to our code:
>
> MeshCommunication().gather_neighboring_elements(cast_ref<DistributedMesh
> &>(mesh));
>
> but it seems to take forever to execute that in dbg mode (Even with very
> few grid points i.e. ~100 [8 procs]).
>
> 2)
>
> The mesh looks ok if I output it even without the above line added but I
> trip different asserts in dbg mode depending on the type of element I use
> (I have both implemented).
>
> If I use Tri3 I trip:
>
> Assertion `oldvar == static_cast<Told>(static_cast<Tnew>(oldvar))' failed.
> oldvar = 18446744073709551615
> static_cast<Told>(static_cast<Tnew>(oldvar)) = 4294967295
>
> If I use Quad4 I trip:
>
> Assertion `i < _points.size()' failed.
> i = 0
> _points.size() = 0
>
> The messages are identical on each process.
>
>
> Best
> Barna & Maria
>
> ------------------------------------------------------------
> ------------------
> Check out the vibrant tech community on one of the world's most
> engaging tech sites, SlashDot.org! http://sdm.link/slashdot
> _______________________________________________
> Libmesh-devel mailing list
> Lib...@li...
> https://lists.sourceforge.net/lists/listinfo/libmesh-devel
>
>