Re: [Libmesh-users] transient system of second order

[Roy Stogner <roystgnr@ic...>]

Sorry about the delays, but I'm finally finding some time to work
through old emails now.

On Thu, 21 Feb 2013, Manav Bhatia wrote:

> I am now sure if I understand the concept of "elem_fixed_solution".

That's my fault.  The FEMSystem time integration APIs were
"prematurely optimized", readability really suffered, and we're going
to have to change the APIs soon anyway to enable stabilized nonlinear
mass matrices.

At the moment, in mass_residual():

elem_fixed_solution is equal to some fixed (integrator-dependent)
value of u during your timestep

elem_solution is equal to du/dt at the time of evaluation.
   
> [M] = \int_\Omega w^T \du/dt d\Omega + \int_\Omega L(w)^T \tau du/dt d\Omega
>  
> How does the concept of elem_fixed_solution affect this situation with a mass matrix.

Assuming your L(w) is really L(u)(w) (e.g. due to advection velocity
depending on u) then you'll want to use elem_fixed_solution as the
input for that, or you'll want to yell at me to fix the API right now
so you can start writing to something saner.
---
Roy

Re: [Libmesh-users] possible bug in newton_solver ?

[Roy Stogner <roystgnr@ic...>]

On Thu, 28 Feb 2013, Manav Bhatia wrote:

>  I wanted to run this by the developers to see if this might be an
> unintended behavior:
>
>  Lines 419 in newton_solver.cpp sets the steplength to 1.0, line 420 makes
> a call to line_search which returns the value of step length. However,
> there is no lvalue speficied to the call to line search, so that steplength
> remain 1.0.
>
>  I think the call to line_search should be
>
>  Real steplength = this->line_search(.....);

This is embarrassing.  As best as I can tell from "svn blame", I
committed that bug (a change which ironically enough was being used
for debugging purposes) along with the change set I'd been debugging.

It looks like we are properly modifying the newton_iterate in
line_search itself, though, and the later norm_delta multiplications
and divisions by step_length cancel out, so this is only ruining the
norm_delta printing in verbose mode?

I'll revert it immediately; thank you for catching this!
---
Roy

Re: [Libmesh-users] transient system of second order

[Roy Stogner <roystgnr@ic...>]

On Thu, 21 Feb 2013, Manav Bhatia wrote:

> What is throwing me off is that the class documentation of Euler2_Solver says: 
> 
>  * Euler solves u' = f(theta*u_new + (1-theta)*u_old),
>  * Euler2 solves u' = theta*f(u_new) + (1-theta)*f(u_old)
> 
> which seems to imply that no mass matrix is accounted for.

This is a flaw in the documentation; I'll fix it ASAP.

> Also, while going through the code, I came across the following in euler_solver.C (lines 124-129)
> 
>   // We do a trick here to avoid using a non-1
>   // elem_solution_derivative:
>   context.elem_jacobian *= -1.0;
>   jacobian_computed = _system.mass_residual(jacobian_computed, context) &&
>     jacobian_computed;
>   context.elem_jacobian *= -1.0;
> 
> 
> What is the need to multiply the elem_jacobian by -1.0 both before
> and after the function call? Wouldn't the call to mass_residual
> overwrite the elem_jacobian anyways? 

The call to mass_residual (or anything else) is supposed to *augment*
the Jacobian, rather than writing to a separate jacobian that would
then just be added to the original.  This was my bright idea for
minimizing L1 cache use and is one of the "premature optimization"
examples I was talking about.

A problem with that is, if mass_residual is computing a jacobian
w.r.t. elem_solution, but elem_solution is storing something whose
derivative w.r.t. the real solution isn't 1.0, then the user has to
remember to multiply those jacobian terms by elem_solution_derivative
to take the chain rule into account.  The way we do scaling with
deltat and 1.0 bits was intended to try and avoid this complication
for most codes.
---
Roy

[Libmesh-users] possible bug in newton_solver ?

[Manav Bhatia <bhatiamanav@gm...>]

Hi,

  I wanted to run this by the developers to see if this might be an
unintended behavior:

  Lines 419 in newton_solver.cpp sets the steplength to 1.0, line 420 makes
a call to line_search which returns the value of step length. However,
there is no lvalue speficied to the call to line search, so that steplength
remain 1.0.

  I think the call to line_search should be

  Real steplength = this->line_search(.....);

   Please comment.

Thanks,
Manav

Re: [Libmesh-users] Error Importing Exodus Mesh

[John Peterson <jwpeterson@gm...>]

On Feb 27, 2013, at 5:24 PM, Carlos Acosta <bdm133@...> wrote:

> Thank you very much. That was actually very helpful. I had to save the mesh as quad4 before importing it to libmesh.
> 
> One thing that I found though is that libmesh is not reading the  number of local nodes "n_local_nodes()=0". there for the number of degrees of freedom is zero.
> I am printing mesh info reading the .e file  vs using MeshTools::Generation::build_square (...) and this is what I get.
> 
> *****EXODUS OUTPUT*******
>  Mesh Information:
>   mesh_dimension()=2
>   spatial_dimension()=3
>   n_nodes()=121
>     n_local_nodes()=0
>   n_elem()=100
>     n_local_elem()=0
>     n_active_elem()=100
>   n_subdomains()=1
>   n_partitions()=1
>   n_processors()=1
>   n_threads()=1
>   processor_id()=0
> 
>  EquationSystems
>   n_systems()=1
>    System #0, "Poisson"
>     Type "LinearImplicit"
>     Variables="u" 
>     Finite Element Types="LAGRANGE" 
>     Approximation Orders="FIRST" 
>     n_dofs()=0
>     n_local_dofs()=0
>     n_constrained_dofs()=0
>     n_local_constrained_dofs()=0
>     n_vectors()=1
>     n_matrices()=1
>     DofMap Sparsity
>       Average  On-Processor Bandwidth <= 0
>       Average Off-Processor Bandwidth <= 0
>       Maximum  On-Processor Bandwidth <= 0
>       Maximum Off-Processor Bandwidth <= 0
>     DofMap Constraints
>       Number of DoF Constraints = 0
> *******************
> 
> vs.
> 
> ****************mesh generation output******
> Mesh Information:
>   mesh_dimension()=2
>   spatial_dimension()=3
>   n_nodes()=256
>     n_local_nodes()=256
>   n_elem()=225
>     n_local_elem()=225
>     n_active_elem()=225
>   n_subdomains()=1
>   n_partitions()=1
>   n_processors()=1
>   n_threads()=1
>   processor_id()=0
> 
>  EquationSystems
>   n_systems()=1
>    System #0, "Poisson"
>     Type "LinearImplicit"
>     Variables="u" 
>     Finite Element Types="LAGRANGE" 
>     Approximation Orders="FIRST" 
>     n_dofs()=256
>     n_local_dofs()=256
>     n_constrained_dofs()=0
>     n_local_constrained_dofs()=0
>     n_vectors()=1
>     n_matrices()=1
>     DofMap Sparsity
>       Average  On-Processor Bandwidth <= 8.26562
>       Average Off-Processor Bandwidth <= 0
>       Maximum  On-Processor Bandwidth <= 9
>       Maximum Off-Processor Bandwidth <= 0
>     DofMap Constraints
>       Number of DoF Constraints = 0
> ************************************************
> 
> Is there an extra step I need to do to read the exodus mesh???

Call mesh.prepare_for_use() maybe?


> 
> Thank you in advance.
> 
> 
>                                         
> 
> On Wed, Feb 27, 2013 at 1:06 PM, John Peterson <jwpeterson@...> wrote:
>> 
>> 
>> 
>> 
>> On Feb 27, 2013, at 12:23 PM, "Paul T. Bauman" <ptbauman@...> wrote:
>> 
>> > On Wed, Feb 27, 2013 at 11:51 AM, Carlos Acosta <bdm133@...> wrote:
>> >
>> >> but when I import an exodus mesh I get errors
>> >> recognizing the exodus format.
>> >> It says:
>> >>
>> >>
>> >> ERROR! Unrecognized element type_str: SHELL4
>> >
>> > Hmm, this is strange because Exodus reader should understand SHELL4.
>> 
>> But it doesn't quote work right: shell4's actually have six "sides" and so some of our side id mapping stuff doesn't work.  Definitely stick to Quad4's.  I think I will patch this to go back to making shell4 throw an error, maybe something like "please export quad4 from cubit instead"...
>> 
>> 
>> 
>> > Could
>> > you rerun with METHOD=dbg and see if an assert is tripped? How big is your
>> > mesh?
>> 
>> It's possible he's using an old version, maybe a sourceforge SVN version?
>> 
>> 
>> > ------------------------------------------------------------------------------
>> > Everyone hates slow websites. So do we.
>> > Make your web apps faster with AppDynamics
>> > Download AppDynamics Lite for free today:
>> > http://p.sf.net/sfu/appdyn_d2d_feb
>> > _______________________________________________
>> > Libmesh-users mailing list
>> > Libmesh-users@...
>> > https://lists.sourceforge.net/lists/listinfo/libmesh-users
> 
> 
> 
> -- 
> 
> Carlos Acosta Berlinghieri
> Graduate Research Assistant
> Department of Biomedical Engineering Computational Center
> University of Texas at San Antonio
> (210)458-4680 (lab phone) 

Re: [Libmesh-users] Error Importing Exodus Mesh

[Carlos Acosta <bdm133@my...>]

Thank you very much. That was actually very helpful. I had to save the mesh
as quad4 before importing it to libmesh.

One thing that I found though is that libmesh is not reading the  number of
local nodes "n_local_nodes()=0". there for the number of degrees of freedom
is zero.
I am printing mesh info reading the .e file  vs using
MeshTools::Generation::build_square (...) and this is what I get.

*****EXODUS OUTPUT*******
 Mesh Information:
  mesh_dimension()=2
  spatial_dimension()=3
  n_nodes()=121
    n_local_nodes()=0
  n_elem()=100
    n_local_elem()=0
    n_active_elem()=100
  n_subdomains()=1
  n_partitions()=1
  n_processors()=1
  n_threads()=1
  processor_id()=0

 EquationSystems
  n_systems()=1
   System #0, "Poisson"
    Type "LinearImplicit"
    Variables="u"
    Finite Element Types="LAGRANGE"
    Approximation Orders="FIRST"
    n_dofs()=0
    n_local_dofs()=0
    n_constrained_dofs()=0
    n_local_constrained_dofs()=0
    n_vectors()=1
    n_matrices()=1
    DofMap Sparsity
      Average  On-Processor Bandwidth <= 0
      Average Off-Processor Bandwidth <= 0
      Maximum  On-Processor Bandwidth <= 0
      Maximum Off-Processor Bandwidth <= 0
    DofMap Constraints
      Number of DoF Constraints = 0
*******************

vs.

****************mesh generation output******
Mesh Information:
  mesh_dimension()=2
  spatial_dimension()=3
  n_nodes()=256
    n_local_nodes()=256
  n_elem()=225
    n_local_elem()=225
    n_active_elem()=225
  n_subdomains()=1
  n_partitions()=1
  n_processors()=1
  n_threads()=1
  processor_id()=0

 EquationSystems
  n_systems()=1
   System #0, "Poisson"
    Type "LinearImplicit"
    Variables="u"
    Finite Element Types="LAGRANGE"
    Approximation Orders="FIRST"
    n_dofs()=256
    n_local_dofs()=256
    n_constrained_dofs()=0
    n_local_constrained_dofs()=0
    n_vectors()=1
    n_matrices()=1
    DofMap Sparsity
      Average  On-Processor Bandwidth <= 8.26562
      Average Off-Processor Bandwidth <= 0
      Maximum  On-Processor Bandwidth <= 9
      Maximum Off-Processor Bandwidth <= 0
    DofMap Constraints
      Number of DoF Constraints = 0
************************************************

Is there an extra step I need to do to read the exodus mesh???

Thank you in advance.




On Wed, Feb 27, 2013 at 1:06 PM, John Peterson <jwpeterson@...> wrote:

>
>
>
>
> On Feb 27, 2013, at 12:23 PM, "Paul T. Bauman" <ptbauman@...> wrote:
>
> > On Wed, Feb 27, 2013 at 11:51 AM, Carlos Acosta <bdm133@...>
> wrote:
> >
> >> but when I import an exodus mesh I get errors
> >> recognizing the exodus format.
> >> It says:
> >>
> >>
> >> ERROR! Unrecognized element type_str: SHELL4
> >
> > Hmm, this is strange because Exodus reader should understand SHELL4.
>
> But it doesn't quote work right: shell4's actually have six "sides" and so
> some of our side id mapping stuff doesn't work.  Definitely stick to
> Quad4's.  I think I will patch this to go back to making shell4 throw an
> error, maybe something like "please export quad4 from cubit instead"...
>
>
>
> > Could
> > you rerun with METHOD=dbg and see if an assert is tripped? How big is
> your
> > mesh?
>
> It's possible he's using an old version, maybe a sourceforge SVN version?
>
>
> >
> ------------------------------------------------------------------------------
> > Everyone hates slow websites. So do we.
> > Make your web apps faster with AppDynamics
> > Download AppDynamics Lite for free today:
> > http://p.sf.net/sfu/appdyn_d2d_feb
> > _______________________________________________
> > Libmesh-users mailing list
> > Libmesh-users@...
> > https://lists.sourceforge.net/lists/listinfo/libmesh-users
>



-- 

Carlos Acosta Berlinghieri
Graduate Research Assistant
Department of Biomedical Engineering Computational Center
University of Texas at San Antonio
(210)458-4680 (lab phone)

Re: [Libmesh-users] Error Importing Exodus Mesh

[John Peterson <jwpeterson@gm...>]

On Feb 27, 2013, at 12:23 PM, "Paul T. Bauman" <ptbauman@...> wrote:

> On Wed, Feb 27, 2013 at 11:51 AM, Carlos Acosta <bdm133@...> wrote:
> 
>> but when I import an exodus mesh I get errors
>> recognizing the exodus format.
>> It says:
>> 
>> 
>> ERROR! Unrecognized element type_str: SHELL4
> 
> Hmm, this is strange because Exodus reader should understand SHELL4.

But it doesn't quote work right: shell4's actually have six "sides" and so some of our side id mapping stuff doesn't work.  Definitely stick to Quad4's.  I think I will patch this to go back to making shell4 throw an error, maybe something like "please export quad4 from cubit instead"...



> Could
> you rerun with METHOD=dbg and see if an assert is tripped? How big is your
> mesh?

It's possible he's using an old version, maybe a sourceforge SVN version?


> ------------------------------------------------------------------------------
> Everyone hates slow websites. So do we.
> Make your web apps faster with AppDynamics
> Download AppDynamics Lite for free today:
> http://p.sf.net/sfu/appdyn_d2d_feb
> _______________________________________________
> Libmesh-users mailing list
> Libmesh-users@...
> https://lists.sourceforge.net/lists/listinfo/libmesh-users

Re: [Libmesh-users] Error Importing Exodus Mesh

[David Knezevic <dknezevic@se...>]

Is there a difference in Cubit between QUAD4 and SHELL4?



On 02/27/2013 01:50 PM, John Peterson wrote:
> Yes you need to export QUAD4 from cubit instead of SHELL4 which is the default.
>
> --
> John
>
> On Feb 27, 2013, at 11:51 AM, Carlos Acosta <bdm133@...> wrote:
>
>> Hello.
>>
>> I am trying to run one of the first examples. The one solving a Poisson
>> problem.
>> The example runs great when the mesh is generated from Libmesh
>> "MeshTools::Generation::..."  but when I import an exodus mesh I get errors
>> recognizing the exodus format.
>> It says:
>>
>>
>> ERROR! Unrecognized element type_str: SHELL4
>> Stack frames: 6
>> 0: libMesh::print_trace(std::ostream&)
>> 1: libMesh::ExodusII_IO_Helper::ElementMaps::assign_conversion(std::string)
>> 2: libMesh::ExodusII_IO::read(std::string const&)
>> 3: main
>> 4: __libc_start_main
>> 5: ./exemplo_1.2-opt() [0x422049]
>> [0] src/mesh/exodusII_io_helper.C, line 1596, compiled Dec  3 2012 at
>> 19:00:45
>> terminate called after throwing an instance of 'libMesh::LogicError'
>>   what():  Error in libMesh internal logic
>> [cacosta-desktop:16352] *** Process received signal ***
>> [cacosta-desktop:16352] Signal: Aborted (6)
>> [cacosta-desktop:16352] Signal code:  (-6)
>> [cacosta-desktop:16352] [ 0] /lib/libpthread.so.0(+0xf8f0) [0x7f9b9aaa78f0]
>> [cacosta-desktop:16352] [ 1] /lib/libc.so.6(gsignal+0x35) [0x7f9b99d99b25]
>> [cacosta-desktop:16352] [ 2] /lib/libc.so.6(abort+0x180) [0x7f9b99d9d670]
>> [cacosta-desktop:16352] [ 3]
>> /usr/lib/libstdc++.so.6(_ZN9__gnu_cxx27__verbose_terminate_handlerEv+0x115)
>> [0x7f9b9d51d8e5]
>> [cacosta-desktop:16352] [ 4] /usr/lib/libstdc++.so.6(+0xcad16)
>> [0x7f9b9d51bd16]
>> [cacosta-desktop:16352] [ 5] /usr/lib/libstdc++.so.6(+0xcad43)
>> [0x7f9b9d51bd43]
>> [cacosta-desktop:16352] [ 6] /usr/lib/libstdc++.so.6(+0xcae3e)
>> [0x7f9b9d51be3e]
>> [cacosta-desktop:16352] [ 7]
>> /usr/share/libmesh/lib/x86_64-unknown-linux-gnu_opt/libmesh.so(_ZN7libMesh18ExodusII_IO_Helper11ElementMaps17assign_conversionESs+0x4f0)
>> [0x7f9ba2cf1630]
>> [cacosta-desktop:16352] [ 8]
>> /usr/share/libmesh/lib/x86_64-unknown-linux-gnu_opt/libmesh.so(_ZN7libMesh11ExodusII_IO4readERKSs+0x31c)
>> [0x7f9ba2ceca3c]
>> [cacosta-desktop:16352] [ 9] ./exemplo_1.2-opt(main+0x19f) [0x4269af]
>> [cacosta-desktop:16352] [10] /lib/libc.so.6(__libc_start_main+0xfd)
>> [0x7f9b99d84c4d]
>> [cacosta-desktop:16352] [11] ./exemplo_1.2-opt() [0x422049]
>> [cacosta-desktop:16352] *** End of error message ***
>> Aborted
>>
>> *****
>> Any clues about what might be missing.
>>
>>
>> -- 
>>
>> Carlos Acosta Berlinghieri
>> Graduate Research Assistant
>> Department of Biomedical Engineering Computational Center
>> University of Texas at San Antonio
>> (210)458-4680 (lab phone)
>> ------------------------------------------------------------------------------
>> Everyone hates slow websites. So do we.
>> Make your web apps faster with AppDynamics
>> Download AppDynamics Lite for free today:
>> http://p.sf.net/sfu/appdyn_d2d_feb
>> _______________________________________________
>> Libmesh-users mailing list
>> Libmesh-users@...
>> https://lists.sourceforge.net/lists/listinfo/libmesh-users
> ------------------------------------------------------------------------------
> Everyone hates slow websites. So do we.
> Make your web apps faster with AppDynamics
> Download AppDynamics Lite for free today:
> http://p.sf.net/sfu/appdyn_d2d_feb
> _______________________________________________
> Libmesh-users mailing list
> Libmesh-users@...
> https://lists.sourceforge.net/lists/listinfo/libmesh-users

Re: [Libmesh-users] Error Importing Exodus Mesh

[Paul T. Bauman <ptbauman@gm...>]

On Wed, Feb 27, 2013 at 11:51 AM, Carlos Acosta <bdm133@...> wrote:

> but when I import an exodus mesh I get errors
> recognizing the exodus format.
> It says:
>
>
> ERROR! Unrecognized element type_str: SHELL4
>

Hmm, this is strange because Exodus reader should understand SHELL4. Could
you rerun with METHOD=dbg and see if an assert is tripped? How big is your
mesh?

Re: [Libmesh-users] Error Importing Exodus Mesh

[John Peterson <jwpeterson@gm...>]

Yes you need to export QUAD4 from cubit instead of SHELL4 which is the default.

--
John

On Feb 27, 2013, at 11:51 AM, Carlos Acosta <bdm133@...> wrote:

> Hello.
> 
> I am trying to run one of the first examples. The one solving a Poisson
> problem.
> The example runs great when the mesh is generated from Libmesh
> "MeshTools::Generation::..."  but when I import an exodus mesh I get errors
> recognizing the exodus format.
> It says:
> 
> 
> ERROR! Unrecognized element type_str: SHELL4
> Stack frames: 6
> 0: libMesh::print_trace(std::ostream&)
> 1: libMesh::ExodusII_IO_Helper::ElementMaps::assign_conversion(std::string)
> 2: libMesh::ExodusII_IO::read(std::string const&)
> 3: main
> 4: __libc_start_main
> 5: ./exemplo_1.2-opt() [0x422049]
> [0] src/mesh/exodusII_io_helper.C, line 1596, compiled Dec  3 2012 at
> 19:00:45
> terminate called after throwing an instance of 'libMesh::LogicError'
>  what():  Error in libMesh internal logic
> [cacosta-desktop:16352] *** Process received signal ***
> [cacosta-desktop:16352] Signal: Aborted (6)
> [cacosta-desktop:16352] Signal code:  (-6)
> [cacosta-desktop:16352] [ 0] /lib/libpthread.so.0(+0xf8f0) [0x7f9b9aaa78f0]
> [cacosta-desktop:16352] [ 1] /lib/libc.so.6(gsignal+0x35) [0x7f9b99d99b25]
> [cacosta-desktop:16352] [ 2] /lib/libc.so.6(abort+0x180) [0x7f9b99d9d670]
> [cacosta-desktop:16352] [ 3]
> /usr/lib/libstdc++.so.6(_ZN9__gnu_cxx27__verbose_terminate_handlerEv+0x115)
> [0x7f9b9d51d8e5]
> [cacosta-desktop:16352] [ 4] /usr/lib/libstdc++.so.6(+0xcad16)
> [0x7f9b9d51bd16]
> [cacosta-desktop:16352] [ 5] /usr/lib/libstdc++.so.6(+0xcad43)
> [0x7f9b9d51bd43]
> [cacosta-desktop:16352] [ 6] /usr/lib/libstdc++.so.6(+0xcae3e)
> [0x7f9b9d51be3e]
> [cacosta-desktop:16352] [ 7]
> /usr/share/libmesh/lib/x86_64-unknown-linux-gnu_opt/libmesh.so(_ZN7libMesh18ExodusII_IO_Helper11ElementMaps17assign_conversionESs+0x4f0)
> [0x7f9ba2cf1630]
> [cacosta-desktop:16352] [ 8]
> /usr/share/libmesh/lib/x86_64-unknown-linux-gnu_opt/libmesh.so(_ZN7libMesh11ExodusII_IO4readERKSs+0x31c)
> [0x7f9ba2ceca3c]
> [cacosta-desktop:16352] [ 9] ./exemplo_1.2-opt(main+0x19f) [0x4269af]
> [cacosta-desktop:16352] [10] /lib/libc.so.6(__libc_start_main+0xfd)
> [0x7f9b99d84c4d]
> [cacosta-desktop:16352] [11] ./exemplo_1.2-opt() [0x422049]
> [cacosta-desktop:16352] *** End of error message ***
> Aborted
> 
> *****
> Any clues about what might be missing.
> 
> 
> -- 
> 
> Carlos Acosta Berlinghieri
> Graduate Research Assistant
> Department of Biomedical Engineering Computational Center
> University of Texas at San Antonio
> (210)458-4680 (lab phone)
> ------------------------------------------------------------------------------
> Everyone hates slow websites. So do we.
> Make your web apps faster with AppDynamics
> Download AppDynamics Lite for free today:
> http://p.sf.net/sfu/appdyn_d2d_feb
> _______________________________________________
> Libmesh-users mailing list
> Libmesh-users@...
> https://lists.sourceforge.net/lists/listinfo/libmesh-users

[Libmesh-users] Error Importing Exodus Mesh

[Carlos Acosta <bdm133@my...>]

Hello.

I am trying to run one of the first examples. The one solving a Poisson
problem.
The example runs great when the mesh is generated from Libmesh
"MeshTools::Generation::..."  but when I import an exodus mesh I get errors
recognizing the exodus format.
It says:


ERROR! Unrecognized element type_str: SHELL4
Stack frames: 6
0: libMesh::print_trace(std::ostream&)
1: libMesh::ExodusII_IO_Helper::ElementMaps::assign_conversion(std::string)
2: libMesh::ExodusII_IO::read(std::string const&)
3: main
4: __libc_start_main
5: ./exemplo_1.2-opt() [0x422049]
[0] src/mesh/exodusII_io_helper.C, line 1596, compiled Dec  3 2012 at
19:00:45
terminate called after throwing an instance of 'libMesh::LogicError'
  what():  Error in libMesh internal logic
[cacosta-desktop:16352] *** Process received signal ***
[cacosta-desktop:16352] Signal: Aborted (6)
[cacosta-desktop:16352] Signal code:  (-6)
[cacosta-desktop:16352] [ 0] /lib/libpthread.so.0(+0xf8f0) [0x7f9b9aaa78f0]
[cacosta-desktop:16352] [ 1] /lib/libc.so.6(gsignal+0x35) [0x7f9b99d99b25]
[cacosta-desktop:16352] [ 2] /lib/libc.so.6(abort+0x180) [0x7f9b99d9d670]
[cacosta-desktop:16352] [ 3]
/usr/lib/libstdc++.so.6(_ZN9__gnu_cxx27__verbose_terminate_handlerEv+0x115)
[0x7f9b9d51d8e5]
[cacosta-desktop:16352] [ 4] /usr/lib/libstdc++.so.6(+0xcad16)
[0x7f9b9d51bd16]
[cacosta-desktop:16352] [ 5] /usr/lib/libstdc++.so.6(+0xcad43)
[0x7f9b9d51bd43]
[cacosta-desktop:16352] [ 6] /usr/lib/libstdc++.so.6(+0xcae3e)
[0x7f9b9d51be3e]
[cacosta-desktop:16352] [ 7]
/usr/share/libmesh/lib/x86_64-unknown-linux-gnu_opt/libmesh.so(_ZN7libMesh18ExodusII_IO_Helper11ElementMaps17assign_conversionESs+0x4f0)
[0x7f9ba2cf1630]
[cacosta-desktop:16352] [ 8]
/usr/share/libmesh/lib/x86_64-unknown-linux-gnu_opt/libmesh.so(_ZN7libMesh11ExodusII_IO4readERKSs+0x31c)
[0x7f9ba2ceca3c]
[cacosta-desktop:16352] [ 9] ./exemplo_1.2-opt(main+0x19f) [0x4269af]
[cacosta-desktop:16352] [10] /lib/libc.so.6(__libc_start_main+0xfd)
[0x7f9b99d84c4d]
[cacosta-desktop:16352] [11] ./exemplo_1.2-opt() [0x422049]
[cacosta-desktop:16352] *** End of error message ***
Aborted

*****
Any clues about what might be missing.


-- 

Carlos Acosta Berlinghieri
Graduate Research Assistant
Department of Biomedical Engineering Computational Center
University of Texas at San Antonio
(210)458-4680 (lab phone)

[Libmesh-users] Fwd: Compiling problems

[Cody Permann <codypermann@gm...>]

I have a user trying to configure libMesh on a Linux cluster but it keeps
running into problems with TBB.  He's not alone as I have also had a few
other reports of people who still receive error even after adding the
--diable-tbb flag to the configure line.  Are we being overly-agressive
when snooping for tbb?  I haven't really dived into this much yet but just
wanted to see if anyone is aware of a quick fix.

The compile errros are attached at the bottom of this thread.

Thanks,
Cody


---------- Forwarded message ----------
From: Capps, Nathan Allen <ncapps@...>
Date: Wed, Feb 27, 2013 at 8:30 AM
Subject: RE: Compiling problems
To: "Permann, Cody J" <cody.permann@...>
Cc: "marie.backman@..." <marie.backman@...>


 I disabled tbb as you suggested.  It seem to produce a similar error.

 Nathan
 ------------------------------
*From:* Permann, Cody J [cody.permann@...]
*Sent:* Wednesday, February 27, 2013 10:17 AM

*To:* Capps, Nathan Allen
*Subject:* Re: Compiling problems


If you take a look at the "build_libmesh_moose.sh" script.  It essentially
calls libMesh's configure with a couple of switches we need to build moose.
 You can simply add more flags to that script, in this case "--disable-tbb".

 You can always run "./configure --help" in libmesh to see what flags and
options are available.  It's possible you might have to turn off cppthreads
too depending on how your system is setup.

 Cody


On Tue, Feb 26, 2013 at 5:29 PM, Capps, Nathan Allen <ncapps@...> wrote:

>  Thank you for your input and how do I turn off packages which are not
> needed and how do I turn off Tbb
>
>  Nathan
>
> Sent from my iPhone
>
> On Feb 26, 2013, at 6:28 PM, "Permann, Cody J" <cody.permann@...>
> wrote:
>
>   Well - all of your problems are pointing to an issue with TBB (Intel's
> Threaded Building Blocks).  It looks like libMesh is detecting it but
> something isn't quite right with that package.  You might want to just
> disable (--disable-tbb).  Yeah - you won't be able to run threaded but
> that's probably not much of a concern until you are running jobs on several
> thousand processors.
>
>  In general, there are a lot of packages that libMesh looks for and
> configures if they are on the system.  You don't need them all so you can
> either try to figure out why it's not working or you can just turn it off.
>
>  BTW - if you want to send output (which is useful) just send it as plain
> text :)
>
>  Cody
>
>
> On Tue, Feb 26, 2013 at 4:10 PM, Capps, Nathan Allen <ncapps@...>wrote:
>
>>  Cody,
>>
>>  I apologize for this being a week later. Here is the output when
>> compiling lib mesh. I sent you the whole output in hopes to essentially
>> help you help me lol.  Also I have cc a Brian's new post doc. She is going
>> to be using the MOOSE system and is helping me get MOOSE compiled on are
>> local cluster.  If there is information you need, just let me know.
>>
>>  Nathan
>>   ------------------------------
>> *From:* Permann, Cody J [cody.permann@...]
>> *Sent:* Wednesday, February 20, 2013 3:58 PM
>> *To:* Capps, Nathan Allen
>> *Subject:* Re: Compiling problems
>>
>>    Are you using our libmesh build script?  Is there more to this
>> message? It tells me which package is broken but not why.
>>
>>  Cody
>>
>>
>> On Wed, Feb 20, 2013 at 1:35 PM, Capps, Nathan Allen <ncapps@...>wrote:
>>
>>>  Cody,
>>>
>>>  I have been trying to compile the Moose system on my local cluster,
>>> and I have been running into a few errors in libmesh listed below.  Also I
>>> have the petcs-3.3-p4. Do you have any clues into what the problem might be?
>>>
>>>  collect2: ld returned 1 exit status
>>> make[1]: *** [getpot_parse-opt] Error 1
>>> make[1]: Leaving directory `/home/ncapps/projects/trunk/libmesh/build'
>>> make: *** [install-recursive] Error 1
>>>
>>>  Nathan
>>>
>>
>>
>

Re: [Libmesh-users] scalar variables

[John Peterson <jwpeterson@gm...>]

I think using higher order scalar variables is more or less equivalent to using more than one first-order scalar.  

It can be a convenient way of grouping together scalar variables which are logically related in some way.

--
John

On Feb 26, 2013, at 8:53 AM, Manav Bhatia <bhatiamanav@...> wrote:

> Hi,
> 
>   I am trying to understand the SCALAR variable/fe types.
> 
>   The FE class of SCALAR seems to return 0 for the shape functions in all
> cases.  At the same time, the project solution method/class uses
> n_components of a SCALAR variable, which in turn is the order of the scalar
> variable.
> 
>   So, if scalar variable can have higher orders with constant (0.0) shape
> function, in which situations can n_components be useful as opposed to
> adding multiple scalar variables to a system?
> 
>   Any comments would be apprecited.
> 
> Thanks,
> Manav
> ------------------------------------------------------------------------------
> Everyone hates slow websites. So do we.
> Make your web apps faster with AppDynamics
> Download AppDynamics Lite for free today:
> http://p.sf.net/sfu/appdyn_d2d_feb
> _______________________________________________
> Libmesh-users mailing list
> Libmesh-users@...
> https://lists.sourceforge.net/lists/listinfo/libmesh-users

[Libmesh-users] scalar variables

[Manav Bhatia <bhatiamanav@gm...>]

Hi,

   I am trying to understand the SCALAR variable/fe types.

   The FE class of SCALAR seems to return 0 for the shape functions in all
cases.  At the same time, the project solution method/class uses
n_components of a SCALAR variable, which in turn is the order of the scalar
variable.

   So, if scalar variable can have higher orders with constant (0.0) shape
function, in which situations can n_components be useful as opposed to
adding multiple scalar variables to a system?

   Any comments would be apprecited.

Thanks,
Manav

[Libmesh-users] Limited timestep size in MemorySolutionHistory

[Robert <libmesh@ro...>]

Hello,

after looking into it, I think that the MemorySolutionHistory is only
usable for timesteps greater than TOLERANCE (which is 1e-6 here):
> https://github.com/libMesh/libmesh/blob/master/src/solvers/memory_solution_history.C#L108

Is it possible to modify this class to allow smaller timesteps?

Regards,
Robert

Re: [Libmesh-users] transient nonlinear implicit system

[Subramanya Gautam Sadasiva <ssadasiv@pu...>]

Hi, I was looking through the code and I realized that a clear method is called for the nonlinear_solver irrespective of whether it is called through the a NonlinearImplicitSystem or a TransientNonlinearImplicitSystem. This strikes me as a bit wasteful. Is there are a reason for releasing all the memory allocated to the solver, and reconstructing it? I can override it but it seems to be a quite a tedious way to go about it. Is there an easy workaround for this? I am quite new to c++ so it doesn't strike me as a very easy thing to do. 
Thanks, 
Subramanya Sadasiva 

"But memory eventually fades. Turbulences damp out, internal strains yield to plastic flow, concentration inhomogeneities diffuse to uniformity. Systems tend to subside to very simple states,independent of their specific history."
Herbert Callen . Thermodynamics and an Introduction to Thermostatics.

----- Original Message -----
From: "Subramanya Gautam Sadasiva" <ssadasiv@...>
To: "libmesh-users" <libmesh-users@...>
Sent: Monday, February 25, 2013 3:39:27 PM
Subject: [Libmesh-users] transient nonlinear implicit system

Hi, 
I am trying to use SNESSetPrefix to set a name to the solver to set options for the solver. I do this by deriving from TransientNonlinearImplicitSystem and I do 
 
    PetscDMNonlinearSolver<Number>*  petsc_solver = 
    new PetscDMNonlinearSolver<Number> (*this); 
  
  // Attach the R & J calculation object
  petsc_solver->residual_object = this;
  petsc_solver->jacobian_object = this;
  // Attach the bounds calculation object 
  petsc_solver->bounds_object = this;
  
  // Get a petsc error . this is an alias for a integer 
  PetscErrorCode ierr; 
  
  petsc_solver->PetscNonlinearSolver<Number>::init(); 
 
  // The prefix is ch_solve_. I need to figure out how to 
  // put it into the documentation 
  ierr = SNESSetOptionsPrefix(petsc_solver->_snes,"ch_solve_");
 

in the constructor. This setting lasts for the first call to solve(). However - when I call the solver the second time, all the options are lost. Am I doing something wrong? Do I need to do something else to prevent the solver from being completely re-initialized. 
Thanks, 

Subramanya Sadasiva 

"But memory eventually fades. Turbulences damp out, internal strains yield to plastic flow, concentration inhomogeneities diffuse to uniformity. Systems tend to subside to very simple states,independent of their specific history."
Herbert Callen . Thermodynamics and an Introduction to Thermostatics.

------------------------------------------------------------------------------
Everyone hates slow websites. So do we.
Make your web apps faster with AppDynamics
Download AppDynamics Lite for free today:
http://p.sf.net/sfu/appdyn_d2d_feb
_______________________________________________
Libmesh-users mailing list
Libmesh-users@...
https://lists.sourceforge.net/lists/listinfo/libmesh-users

[Libmesh-users] transient nonlinear implicit system

[Subramanya Gautam Sadasiva <ssadasiv@pu...>]

Hi, 
I am trying to use SNESSetPrefix to set a name to the solver to set options for the solver. I do this by deriving from TransientNonlinearImplicitSystem and I do 
 
    PetscDMNonlinearSolver<Number>*  petsc_solver = 
    new PetscDMNonlinearSolver<Number> (*this); 
  
  // Attach the R & J calculation object
  petsc_solver->residual_object = this;
  petsc_solver->jacobian_object = this;
  // Attach the bounds calculation object 
  petsc_solver->bounds_object = this;
  
  // Get a petsc error . this is an alias for a integer 
  PetscErrorCode ierr; 
  
  petsc_solver->PetscNonlinearSolver<Number>::init(); 
 
  // The prefix is ch_solve_. I need to figure out how to 
  // put it into the documentation 
  ierr = SNESSetOptionsPrefix(petsc_solver->_snes,"ch_solve_");
 

in the constructor. This setting lasts for the first call to solve(). However - when I call the solver the second time, all the options are lost. Am I doing something wrong? Do I need to do something else to prevent the solver from being completely re-initialized. 
Thanks, 

Subramanya Sadasiva 

"But memory eventually fades. Turbulences damp out, internal strains yield to plastic flow, concentration inhomogeneities diffuse to uniformity. Systems tend to subside to very simple states,independent of their specific history."
Herbert Callen . Thermodynamics and an Introduction to Thermostatics.

[Libmesh-users] getting u and u_dot from transient solver

[Manav Bhatia <bhatiamanav@gm...>]

Hi, 

   I am attempting to implement my SUPG solver for Euler equations within libMesh via the FEMSystem API. A couple of questions: 

--  I am assuming that when elem_time_derivative is called, elem_solution in the FEMContext stores the u. Similarly, when the mass_residual is called, elem_solution stores the u_dot and elem_fixed_solution stores the u. Is that correct? 

-- How do I get the element u_dot vector when elem_time_derivative is called? 

    I would appreciate any help. 

Thanks,
Manav

Re: [Libmesh-users] Forcing the solver to use a particular solver.

[Subramanya Gautam Sadasiva <ssadasiv@pu...>]

Sorry about this mail . I see there is a really stupid error in what I am doing. 

Subramanya Sadasiva 

"But memory eventually fades. Turbulences damp out, internal strains yield to plastic flow, concentration inhomogeneities diffuse to uniformity. Systems tend to subside to very simple states,independent of their specific history."
Herbert Callen . Thermodynamics and an Introduction to Thermostatics.

----- Original Message -----
From: "Subramanya Gautam Sadasiva" <ssadasiv@...>
To: "libmesh-users" <libmesh-users@...>
Sent: Friday, February 22, 2013 2:22:42 PM
Subject: [Libmesh-users] Forcing the solver to use a particular solver.

Hi, 
I want to use the petsc_dm solver without passing --use-petsc-dm at the command line. I did this in the constructor of my class that derives from nonlinearimplicitsystem . 

  AutoPtr<PetscNonlinearSolver<Number> > petsc_solver = 
    AutoPtr<PetscNonlinearSolver<Number> > (new PetscNonlinearSolver<Number> (*this)); 
  
  // Attach the R & J calculation object
  petsc_solver->residual_object = this;
  petsc_solver->jacobian_object = this;
  // Attach the bounds calculation object 
  petsc_solver->bounds_object = this;
 
 // Get a petsc error . this is an alias for a integer 
  PetscErrorCode ierr; 
  
  petsc_solver->PetscNonlinearSolver<Number>::init(); 
 
  // The prefix is ch_solve_. I need to figure out how to 
  // put it into the documentation 
  ierr = SNESSetOptionsPrefix(petsc_solver->_snes,"ch_solve_");
  
  // Set from options file for the solve 
  ierr = SNESSetFromOptions(petsc_solver->_snes); 

  // Now put petsc solver back into nonlinear solver AutoPtr 
  nonlinear_solver = petsc_solver; 
 
This does not seem to work and gives me an error. Is there a better way to do it than to override the nonlinear_solver? 
Thanks,
Subramanya

"But memory eventually fades. Turbulences damp out, internal strains yield to plastic flow, concentration inhomogeneities diffuse to uniformity. Systems tend to subside to very simple states,independent of their specific history."
Herbert Callen . Thermodynamics and an Introduction to Thermostatics.


------------------------------------------------------------------------------
Everyone hates slow websites. So do we.
Make your web apps faster with AppDynamics
Download AppDynamics Lite for free today:
http://p.sf.net/sfu/appdyn_d2d_feb
_______________________________________________
Libmesh-users mailing list
Libmesh-users@...
https://lists.sourceforge.net/lists/listinfo/libmesh-users

[Libmesh-users] Forcing the solver to use a particular solver.

[Subramanya Gautam Sadasiva <ssadasiv@pu...>]

Hi, 
I want to use the petsc_dm solver without passing --use-petsc-dm at the command line. I did this in the constructor of my class that derives from nonlinearimplicitsystem . 

  AutoPtr<PetscNonlinearSolver<Number> > petsc_solver = 
    AutoPtr<PetscNonlinearSolver<Number> > (new PetscNonlinearSolver<Number> (*this)); 
  
  // Attach the R & J calculation object
  petsc_solver->residual_object = this;
  petsc_solver->jacobian_object = this;
  // Attach the bounds calculation object 
  petsc_solver->bounds_object = this;
 
 // Get a petsc error . this is an alias for a integer 
  PetscErrorCode ierr; 
  
  petsc_solver->PetscNonlinearSolver<Number>::init(); 
 
  // The prefix is ch_solve_. I need to figure out how to 
  // put it into the documentation 
  ierr = SNESSetOptionsPrefix(petsc_solver->_snes,"ch_solve_");
  
  // Set from options file for the solve 
  ierr = SNESSetFromOptions(petsc_solver->_snes); 

  // Now put petsc solver back into nonlinear solver AutoPtr 
  nonlinear_solver = petsc_solver; 
 
This does not seem to work and gives me an error. Is there a better way to do it than to override the nonlinear_solver? 
Thanks,
Subramanya

"But memory eventually fades. Turbulences damp out, internal strains yield to plastic flow, concentration inhomogeneities diffuse to uniformity. Systems tend to subside to very simple states,independent of their specific history."
Herbert Callen . Thermodynamics and an Introduction to Thermostatics.

Re: [Libmesh-users] transient system of second order

[Manav Bhatia <bhatiamanav@gm...>]

	
On Feb 21, 2013, at 3:25 PM, Roy Stogner <roystgnr@...> wrote:

>> 
>>  Additionally, within the existing framework, how does one account for
>> *nonlinear* systems that have a mass matrix
>> 
>>  [M] du/dt = f(u)
> 
> Right now, we allow you to request an "elem_fixed_solution" that can
> be used for either low-order time integration of nonlinear mass
> matrices or for nonlinear stabilization in test functions.  We're
> actually planning a (backwards-incompatible, unfortunately) change in
> this API, to make it less unnecessarily confusing and to allow for
> higher-order time integration of nonlinear mass matrices and/or for
> nonlinear mass matrices with nonlinear stabilization in test
> functions.  If you'd like, I could probably push that change up higher
> on my todo list.


Upon further study, it seems like the elem_fixed_solution stores the element solution in situations where the elem_solution may contain the velocity or delta solution. Is this correct? 

Manav

Re: [Libmesh-users] transient system of second order

[Manav Bhatia <bhatiamanav@gm...>]

On Feb 21, 2013, at 3:25 PM, Roy Stogner <roystgnr@...> wrote:

> 
>>  Additionally, within the existing framework, how does one account for
>> *nonlinear* systems that have a mass matrix
>> 
>>  [M] du/dt = f(u)
> 
> Right now, we allow you to request an "elem_fixed_solution" that can
> be used for either low-order time integration of nonlinear mass
> matrices or for nonlinear stabilization in test functions.  We're
> actually planning a (backwards-incompatible, unfortunately) change in
> this API, to make it less unnecessarily confusing and to allow for
> higher-order time integration of nonlinear mass matrices and/or for
> nonlinear mass matrices with nonlinear stabilization in test
> functions.  If you'd like, I could probably push that change up higher
> on my todo list.
> 

OK, I think I may have started to understand this a little better now (although still a little confused). Please correct me as you see fit.  

I see that the method mass_residual is used in the classes Euler_Solver and Euler2_Solver. 

What is throwing me off is that the class documentation of Euler2_Solver says: 

 * Euler solves u' = f(theta*u_new + (1-theta)*u_old),
 * Euler2 solves u' = theta*f(u_new) + (1-theta)*f(u_old)

which seems to imply that no mass matrix is accounted for. Are the u and f(u) here the solution after discretization, or before discretization? 

My current guess is that the u' in these equations should actually be replaced with the product [M] u' , like one would encounter post-discretization. What is your assessment? 

Also, while going through the code, I came across the following in euler_solver.C (lines 124-129)

  // We do a trick here to avoid using a non-1
  // elem_solution_derivative:
  context.elem_jacobian *= -1.0;
  jacobian_computed = _system.mass_residual(jacobian_computed, context) &&
    jacobian_computed;
  context.elem_jacobian *= -1.0;


What is the need to multiply the elem_jacobian by -1.0 both before and after the function call? Wouldn't the call to mass_residual overwrite the elem_jacobian anyways? 


Thanks,
Manav

Re: [Libmesh-users] Interpolation order

[Manav Bhatia <bhatiamanav@gm...>]

I had looked at DirichletBoundary class in the examples, but I had not
realized that it works for all interpolation basis.

What you describe is very impressive, and I will certaintly look deeper
into the class.

Thanks!
Manav



On Thu, Feb 21, 2013 at 4:17 PM, Roy Stogner <roystgnr@...>wrote:

>
> On Thu, 21 Feb 2013, Manav Bhatia wrote:
>
>     I have a preference for Lagrange interpolation for the structural
>> (and some thermal) applications that I am working on, primarily due
>> to easy application of Dirichlet boundary conditions.
>>
>
> Have you noticed the relatively-new DirichletBoundary class we have?
> Hand that a functor, and it'll do the projections onto non-Lagrange
> spaces for you.  Right now there's still work and testing to be done
> for transient Dirichlet conditions (that functor may get evaluated at
> time t in places where it should be evaluated at t +/- deltat...) but
> for steady Dirichlet conditions it's pretty solid.
>
>
> Perhaps hierarchic polynomials would also offer a similar feature
>> (please correct if I am wrong) due to the presence of node and edge
>> functions whose dof values can be specified, and bubble functions
>> that are zero on the boundaries.
>>
>
> Unless your Dirichlet values are piecewise-linear, they'll have
> non-zero dofs corresponding to those bubble functions.  What our
> DirichletBoundary stuff gives you is interpolation at nodes, then
> constrained L2 projection for edge-interior values, then
> further-constrained L2 projection for face-interior values.
>
>
>    I am curious about your comment "mapped the Lagrange bases on to
>> nodes topologically". Are you implying a software fix of telling the
>> dof that it is associated with a node, or some sort of a
>> mathematical mapping process?
>>
>
> Software.  Suppose you have a quad on the unit square.  Isoparametric
> cubic Lagrange elements would have nodes at (1/3,0) and (2/3,0), each
> with a single degree of freedom.  My suggestion would be to do what we
> do with hierarchics: keep the single node at (1/2,0), give it two
> degrees of freedom, then use some canonical ordering based on node
> positions to determine which DoF is which.  This breaks code which
> assumes that solution(node->dof_number(0,0,**0)) == u(Point(node)), but
> mathematically it can still be the same basis being used.
> ---
> Roy

Re: [Libmesh-users] transient system of second order

[Manav Bhatia <bhatiamanav@gm...>]

Hi Roy,

   In the FEM code that I have written, I have implemented Newmark Method
for first and second order nonlinear systems, with and without mass and
damping matrices. While converting second-order to first-order
systems I create a new sparsity patter data structure where the
sub-matrices inherit their patterns from the original system, or get the
additional nonzeros that show up from the Identity matrices.

   I can look at porting this over to libMesh API, and your idea of
dof_coupling may serve the necessary purpose here. I can look into this.

   I am now sure if I understand the concept of "elem_fixed_solution". For
the Euler equations

  L (u) = 0
  L = d/dt + A_i d/d x_i

 that I have implemented, I get the following equations

\int_\Omega w^T \du/dt d\Omega -
\int_\Omega dw^T/d x_i A_i u d\Omega +
\int_\Lambda w^T flux(u) \Lambda +
\int_\Omega L(w)^T \tau (d/dt + A_i d/dx_i)u d\Omega = 0

This gives the mass matrix
[M] = \int_\Omega w^T \du/dt d\Omega + \int_\Omega L(w)^T \tau du/dt d\Omega

and the rhs function

f(u) =  - (- \int_\Omega dw^T/d x_i A_i u d\Omega + \int_\Lambda w^T
flux(u) \Lambda +\int_\Omega L(w)^T \tau A_i du/dx_i d\Omega  )


How does the concept of elem_fixed_solution affect this situation with a
mass matrix.

Now, for steady solutions, I have successfully time marched by assumed the
mass matrix to be identity, which brings me back to du/dt = f (u). But, I
am curious about time-accurate simulations, where this assumption may not
be accurate.

I would appreciate your comments.

Thanks,
Manav




On Thu, Feb 21, 2013 at 3:25 PM, Roy Stogner <roystgnr@...>wrote:

>
> On Thu, 21 Feb 2013, Manav Bhatia wrote:
>
>   I would appreciate some help in understanding if the existing System
>> classes can handle PDEs (linear/nonlinear) with second order time
>> derivative.
>>
>
> IIRC We have one existing system (NewmarkSystem?) designed for second
> order time derivatives.  I have no idea what its limitations are or
> whether it would be suitable for you.
>
>
> I am going through the FEMSystem documentation and it seems
>> like these are suitable for a system of the form
>>
>>    du/dt = f(u).
>>
>
> Currently, yes.  There's nothing about the APIs that would prevent you
> from adding a UnsteadySolver subclass that integrates u''=f(u), and
> we'd love to get a patch like that, but we don't have any such
> subclass in the library yet.
>
> Worse: if you wanted to handle even more general second-order time
> derivatives, with mixtures of u'' and u' in the same equation, I think
> we'd need some modification to the APIs before we could handle that
> naturally.
>
>
>   Ofcourse, I can change my second order system from
>>
>>    d^2u/dt^2 = f(u)
>>
>>   to
>>
>>    d/dt {u; \tilde{u}} = [0 I; 0 0] {u; \tilde{u}} + {0; f(u)}
>>
>>   but that changes the dimension of my sparse system by a factor of two,
>>
>
> This is probably the easiest way to get a second-order system up and
> running with FEMSystem.
>
>
> and I am not sure if the FEMSystem takes into account the modification in
>> the sparsity pattern.
>>
>
> You do the add_variable() for your auxilliary u' variable in your
> system initialization, and the sparsity pattern will be modified to
> assume that every residual term corresponding to a u degree of freedom
> on an element is coupled to every u' degree of freedom on that
> element.  This will work; you'll just have unnecessary zeros in your
> sparse matrix or you'll have to use our DofMap::_dof_coupling API to
> avoid them.
>
>
>   Additionally, within the existing framework, how does one account for
>> *nonlinear* systems that have a mass matrix
>>
>>   [M] du/dt = f(u)
>>
>
> Right now, we allow you to request an "elem_fixed_solution" that can
> be used for either low-order time integration of nonlinear mass
> matrices or for nonlinear stabilization in test functions.  We're
> actually planning a (backwards-incompatible, unfortunately) change in
> this API, to make it less unnecessarily confusing and to allow for
> higher-order time integration of nonlinear mass matrices and/or for
> nonlinear mass matrices with nonlinear stabilization in test
> functions.  If you'd like, I could probably push that change up higher
> on my todo list.
>
>
>   Extending this question forward, how about PDEs with second order time
>> derivative?
>>
>
> Now that's a good question.  I've only done second-order time
> integration by breaking it up into first-order systems.  I'd like to
> support direct second-order integrators, but I'd need help designing
> the API so as to make it as broadly applicable as possible.
> ---
> Roy
>

Re: [Libmesh-users] Interpolation order

[Roy Stogner <roystgnr@ic...>]

On Thu, 21 Feb 2013, Manav Bhatia wrote:

>    I have a preference for Lagrange interpolation for the structural
> (and some thermal) applications that I am working on, primarily due
> to easy application of Dirichlet boundary conditions.

Have you noticed the relatively-new DirichletBoundary class we have?
Hand that a functor, and it'll do the projections onto non-Lagrange
spaces for you.  Right now there's still work and testing to be done
for transient Dirichlet conditions (that functor may get evaluated at
time t in places where it should be evaluated at t +/- deltat...) but
for steady Dirichlet conditions it's pretty solid.

> Perhaps hierarchic polynomials would also offer a similar feature
> (please correct if I am wrong) due to the presence of node and edge
> functions whose dof values can be specified, and bubble functions
> that are zero on the boundaries. 

Unless your Dirichlet values are piecewise-linear, they'll have
non-zero dofs corresponding to those bubble functions.  What our
DirichletBoundary stuff gives you is interpolation at nodes, then
constrained L2 projection for edge-interior values, then
further-constrained L2 projection for face-interior values.

>    I am curious about your comment "mapped the Lagrange bases on to
> nodes topologically". Are you implying a software fix of telling the
> dof that it is associated with a node, or some sort of a
> mathematical mapping process?

Software.  Suppose you have a quad on the unit square.  Isoparametric
cubic Lagrange elements would have nodes at (1/3,0) and (2/3,0), each
with a single degree of freedom.  My suggestion would be to do what we
do with hierarchics: keep the single node at (1/2,0), give it two
degrees of freedom, then use some canonical ordering based on node
positions to determine which DoF is which.  This breaks code which
assumes that solution(node->dof_number(0,0,0)) == u(Point(node)), but
mathematically it can still be the same basis being used.
---
Roy