On Mar 16, 2009, at 12:58 PM, Roy Stogner wrote:
> Of course, sometimes we're happy with false negatives. If you're
> looping over all elements, then your technique will catch every
> boundary node eventually, it just won't catch them from each element
> they're in. If you're setting non-penalty Dirichlet conditions, for
> example, that's good enough.
That depends on how you're setting that condition. If you're hardcore
going in and modifying the global matrix then you're correct. If
you're modifying the element matrix _before_ adding it into the global
matrix (as is often done with penalty formulations.... ie Example 3)
then that won't work (because you never modify the element matrix for
that element that has a node but no side on the boundary).
Just thought I would point that out for anyone out there using the SLT
(Standard Libmesh Template... ie modifying an example) and wanting to
set non-penalty conditions....
Derek

On Mon, 16 Mar 2009, Andrea Hawkins wrote:
> I am just wondering if there is a simple way to check if a node is on the
> boundary.
Only if you've assigned a boundary id to that part of the boundary;
then you can just call BoundaryInfo::boundary_id(node) and check to
see if it's not invalid_id.
> The only way I have found I think will work is to first check
> (elem->neighbor(side) == NULL) and then check
> elem.is_node_on_side(i,side).
You'll get false negatives that way. Picture an L-shaped domain
gridded with squares. At that interior corner, there's an element
which has one node but zero sides on the boundary.
Of course, sometimes we're happy with false negatives. If you're
looping over all elements, then your technique will catch every
boundary node eventually, it just won't catch them from each element
they're in. If you're setting non-penalty Dirichlet conditions, for
example, that's good enough.
---
Roy

Hello-
I am just wondering if there is a simple way to check if a node is on the
boundary. The only way I have found I think will work is to first check
(elem->neighbor(side) == NULL) and then check elem.is_node_on_side(i,side).
Did I miss another way?
Thanks!
Andrea

On Sun, 15 Mar 2009, M. Adil Sbai wrote:
> I tried to install LIbMesh under Cygwin in a new laptop (with M$ Windows
> preinstalled) that I got.
>
> Everything went smoothly ... only the linker complained about missing
> definitions !! with the following kind of output:
It's been a long time since I tried libMesh under Cygwin, but it looks
like the linker errors are all confined to two contrib packages.
Would you try configuring with --diable-gzstreams --disable-exodus and
see if that does any better?
---
Roy

If you have MUMPS installed, you can do:
-ksp_type preonly -pc_type lu -mat_type aijmumps
or if you have SuperLu_Dist
-mat_type superlu_dist
both *should* work in parallel with PETSc (unless you compiled MUMPS
with serial only). MUMPS is not as error catching friendly as SuperLU
for zero pivots, etc, but both are reasonably fast sparse-direct solvers.
The key is having a call to the PETSc routine SetMatFromOptions which
will set the MatType from the command line. If you don't have anything,
it just uses what you've set already, or the defaults. This is very
useful for debugging and quickly switching between direct and iterative
solvers. Not sure if libMesh has this call or not, but you can probably
work it in yourself without too much trouble.
FWIW.
Paul
Andrea Hawkins wrote:
> Ah. Ok. Well, I guess I"ll just run it in serial for now. And if I need it
> in parallel, I'll add the option in PetscMatrix and see what happens! =)
>
> Thanks!
> Andrea
>
> On Fri, Mar 13, 2009 at 1:16 PM, John Peterson <jwpeterson@...> wrote:
>
>
>> On Fri, Mar 13, 2009 at 1:15 PM, Andrea Hawkins <andjhawkins@...>
>> wrote:
>>
>>> Well, the suggestion PETSc gives for performing a direct solve is to use
>>>
>> the
>>
>>> options
>>>
>>> -ksp_type preonly -pc_type lu
>>>
>>> But, when I use it I get the error:
>>>
>>> [0]PETSC ERROR: [1]PETSC ERROR: --------------------- Error Message
>>> ------------------------------------
>>> [2]PETSC ERROR: --------------------- Error Message
>>> ------------------------------------
>>> [1]PETSC ERROR: No support for this operation for this object type!
>>> [2]PETSC ERROR: No support for this operation for this object type!
>>> [1]PETSC ERROR: [2]PETSC ERROR: Matrix type mpiaij symbolic LU!
>>> [1]PETSC ERROR: Matrix type mpiaij symbolic LU!
>>> [2]PETSC ERROR:
>>>
>> Sorry, I forgot to mention: it only works in serial.
>>
>> --
>> John
>>
>>
> ------------------------------------------------------------------------------
> Apps built with the Adobe(R) Flex(R) framework and Flex Builder(TM) are
> powering Web 2.0 with engaging, cross-platform capabilities. Quickly and
> easily build your RIAs with Flex Builder, the Eclipse(TM)based development
> software that enables intelligent coding and step-through debugging.
> Download the free 60 day trial. http://p.sf.net/sfu/www-adobe-com
> _______________________________________________
> Libmesh-users mailing list
> Libmesh-users@...
> https://lists.sourceforge.net/lists/listinfo/libmesh-users
>
>

Andrea,
If you want to use LU in parallel, you need to configure petsc with
either MUMPS or SUPERLU_DIST packages. Both these packages will help
you use LU factorization with multiple processors although I am not
exactly sure how they scale (I have not used either of them yet). I
would think the scalability is low since the factorization and
forward/backward solve involve a lot of serial operations that is hard
to parallelize but I could be wrong on this.
By default, Petsc uses Petsc factorization for LU. For more
information, check out MatSolverPackage in Petsc documentation.
Hope that helps.
Vijay
On Fri, Mar 13, 2009 at 1:18 PM, Andrea Hawkins <andjhawkins@...> wrote:
> Ah. Ok. Well, I guess I"ll just run it in serial for now. And if I need it
> in parallel, I'll add the option in PetscMatrix and see what happens! =)
>
> Thanks!
> Andrea
>
> On Fri, Mar 13, 2009 at 1:16 PM, John Peterson <jwpeterson@...> wrote:
>
>> On Fri, Mar 13, 2009 at 1:15 PM, Andrea Hawkins <andjhawkins@...>
>> wrote:
>> > Well, the suggestion PETSc gives for performing a direct solve is to use
>> the
>> > options
>> >
>> > -ksp_type preonly -pc_type lu
>> >
>> > But, when I use it I get the error:
>> >
>> > [0]PETSC ERROR: [1]PETSC ERROR: --------------------- Error Message
>> > ------------------------------------
>> > [2]PETSC ERROR: --------------------- Error Message
>> > ------------------------------------
>> > [1]PETSC ERROR: No support for this operation for this object type!
>> > [2]PETSC ERROR: No support for this operation for this object type!
>> > [1]PETSC ERROR: [2]PETSC ERROR: Matrix type mpiaij symbolic LU!
>> > [1]PETSC ERROR: Matrix type mpiaij symbolic LU!
>> > [2]PETSC ERROR:
>>
>> Sorry, I forgot to mention: it only works in serial.
>>
>> --
>> John
>>
> ------------------------------------------------------------------------------
> Apps built with the Adobe(R) Flex(R) framework and Flex Builder(TM) are
> powering Web 2.0 with engaging, cross-platform capabilities. Quickly and
> easily build your RIAs with Flex Builder, the Eclipse(TM)based development
> software that enables intelligent coding and step-through debugging.
> Download the free 60 day trial. http://p.sf.net/sfu/www-adobe-com
> _______________________________________________
> Libmesh-users mailing list
> Libmesh-users@...
> https://lists.sourceforge.net/lists/listinfo/libmesh-users
>

On Fri, Mar 13, 2009 at 1:13 PM, John Peterson <jwpeterson@...> wrote:
> On Fri, Mar 13, 2009 at 1:09 PM, Andrea Hawkins <andjhawkins@...> wrote:
>> Hello-
>>
>> For debugging purposes, I'm wanting to test out what happens when using a
>> direct solve. However, it appears that for Petsc to do a direct solve the
>> matrix must be a dense matrix. In looking at how libMesh defines
>> PetscMatrix, there does not seem to be an option for this. If I add this
>> option, and use it... about how many things should I expect to crash? Or is
>> there another work around that I just didn't find?
>
> I wouldn't expect too many things to crash ;-)
>
> Would using the LU preconditioner work? Just pass -pc_type lu, and
> any iterative solver should get the solution in 1 iteration.
Oh, and apparently you can also pass
-ksp_type preonly
--
John

On Fri, Mar 13, 2009 at 1:09 PM, Andrea Hawkins <andjhawkins@...> wrote:
> Hello-
>
> For debugging purposes, I'm wanting to test out what happens when using a
> direct solve. However, it appears that for Petsc to do a direct solve the
> matrix must be a dense matrix. In looking at how libMesh defines
> PetscMatrix, there does not seem to be an option for this. If I add this
> option, and use it... about how many things should I expect to crash? Or is
> there another work around that I just didn't find?
I wouldn't expect too many things to crash ;-)
Would using the LU preconditioner work? Just pass -pc_type lu, and
any iterative solver should get the solution in 1 iteration.
--
John

Hello-
For debugging purposes, I'm wanting to test out what happens when using a
direct solve. However, it appears that for Petsc to do a direct solve the
matrix must be a dense matrix. In looking at how libMesh defines
PetscMatrix, there does not seem to be an option for this. If I add this
option, and use it... about how many things should I expect to crash? Or is
there another work around that I just didn't find?
Thanks!
Andrea

On Wed, Mar 11, 2009 at 8:16 AM, Paulo Vieira <paulovieira.ml@...> wrote:
> Hello all!
>
> Here are some more articles that explicitly refer to libmesh for the
> results (found at sciencedirect). You might want to add them to the
> website.
Thanks Paulo!
--
John

On Tue, 10 Mar 2009, Derek Gaston wrote:
> On Mar 10, 2009, at 12:02 PM, Roy Stogner wrote:
>
>> What happens is that phi_face includes all element shape functions,
>> even when they're 0 on the entire face boundary. I guess the choice
>> was between adding a bunch of zeros vs. indexing through an extra
>> indirection layer, and libMesh went with simpler code over
>> negligably better assembly performance.
>
> Also... don't forget non-lagrange elements which could have wacky
> support on boundaries.
Wacky but still well-defined: even for non-Lagrange dofs, if the dof
isn't located on a node on face S, the corresponding shape function
has no support on face S.
But even then you've got to be careful; bases like the monomials may
have no support on interfaces (because they're discontinuous there!)
but you may still want to integrate them there (e.g. to measure that
jump discontinuity). Using every element basis function is just
safer, and even when it's really wasteful (high p elements) it's not
really wasteful (you're just doing it on O(N^2/3) boundary elements,
tops)
> In general I think there's no better way than to just evaluate all
> of the shape functions at the quadrature points and use what comes
> out.
In specific cases there's often a better way, but in general there's
rarely a safer way.
---
Roy

On Mar 10, 2009, at 12:02 PM, Roy Stogner wrote:
> What happens is that phi_face includes all element shape functions,
> even when they're 0 on the entire face boundary. I guess the choice
> was between adding a bunch of zeros vs. indexing through an extra
> indirection layer, and libMesh went with simpler code over
> negligably better assembly performance.
Also... don't forget non-lagrange elements which could have wacky
support on boundaries. In general I think there's no better way than
to just evaluate all of the shape functions at the quadrature points
and use what comes out.
Derek

On Tue, Mar 10, 2009 at 12:04 PM, Andrea Hawkins <andjhawkins@...> wrote:
> Ah. I see. So the extra guys are just set to 0.
>
> Thank you very much!
Yeah, and actually it also works out well for boundary integrals
involving derivatives since (as in your figure above) dphi(0) has
support on the boundary face.
--
John

Ah. I see. So the extra guys are just set to 0.
Thank you very much!
Andrea
On Tue, Mar 10, 2009 at 1:02 PM, Roy Stogner <roystgnr@...>wrote:
>
> On Tue, 10 Mar 2009, Andrea Hawkins wrote:
>
> For discussion I will consider example 5. The loops for test and trial
>> functions run over i,j=0 to phi_face.size() which makes sense, but I'm
>> wondering about the value being added to Ke(i,j) and Fe(i) as I wouldn't
>> think the element dof corresponds to the side-trial function number.
>>
>
> They do correspond.
>
> Is there something that happens in fe_face->reinit() or something?
>>
>
> What happens is that phi_face includes all element shape functions,
> even when they're 0 on the entire face boundary. I guess the choice
> was between adding a bunch of zeros vs. indexing through an extra
> indirection layer, and libMesh went with simpler code over
> negligably better assembly performance.
> ---
> Roy
>

On Tue, 10 Mar 2009, Andrea Hawkins wrote:
> For discussion I will consider example 5. The loops for test and trial
> functions run over i,j=0 to phi_face.size() which makes sense, but I'm
> wondering about the value being added to Ke(i,j) and Fe(i) as I wouldn't
> think the element dof corresponds to the side-trial function number.
They do correspond.
> Is there something that happens in fe_face->reinit() or something?
What happens is that phi_face includes all element shape functions,
even when they're 0 on the entire face boundary. I guess the choice
was between adding a bunch of zeros vs. indexing through an extra
indirection layer, and libMesh went with simpler code over
negligably better assembly performance.
---
Roy

Right... But if you think about a 2D element with local dofs as numbered
here
3 ------------ 2
| |
| |
0 ------------- 1
If the side between dof 2 and 3 is on the boundary, in that loop I don't
think you should be adding to Ke(0,0), right? But you will be. Or have I
completely over thought this?
On Tue, Mar 10, 2009 at 12:54 PM, John Peterson <jwpeterson@...>wrote:
> On Tue, Mar 10, 2009 at 11:51 AM, Andrea Hawkins <andjhawkins@...>
> wrote:
> > Hello-
> >
> > In looking at the examples that include boundary integration, I just had
> a
> > question of how something is working.
> >
> > For discussion I will consider example 5. The loops for test and trial
> > functions run over i,j=0 to phi_face.size() which makes sense, but I'm
> > wondering about the value being added to Ke(i,j) and Fe(i) as I wouldn't
> > think the element dof corresponds to the side-trial function number. Is
> > there something that happens in fe_face->reinit() or something?
>
> All the same basis functions and numbering are used, the quadrature
> rule is changed (I think!)
>
> --
> John
>

On Tue, Mar 10, 2009 at 11:51 AM, Andrea Hawkins <andjhawkins@...> wrote:
> Hello-
>
> In looking at the examples that include boundary integration, I just had a
> question of how something is working.
>
> For discussion I will consider example 5. The loops for test and trial
> functions run over i,j=0 to phi_face.size() which makes sense, but I'm
> wondering about the value being added to Ke(i,j) and Fe(i) as I wouldn't
> think the element dof corresponds to the side-trial function number. Is
> there something that happens in fe_face->reinit() or something?
All the same basis functions and numbering are used, the quadrature
rule is changed (I think!)
--
John

Hello-
In looking at the examples that include boundary integration, I just had a
question of how something is working.
For discussion I will consider example 5. The loops for test and trial
functions run over i,j=0 to phi_face.size() which makes sense, but I'm
wondering about the value being added to Ke(i,j) and Fe(i) as I wouldn't
think the element dof corresponds to the side-trial function number. Is
there something that happens in fe_face->reinit() or something?
Thanks!
Andrea

Hey Roy,
thanks again for your time!
I commented the file writting part out, and changed <Real> type to
<Number>. The compilation was able to get through. Also I collected some
other minor errors due to the complex type. I will list them below
together with the change i made so that the compiler would stop
chattering.
-----------------------------------------------------------------------------
file: src/mesh/ensight_io.C
error:
407 cannot convert ‘std::complex<double>’ to ‘double’
alteration:
379 typedef std::map<int,Real> map_local_soln;
-> typedef std::map<int,Number> map_local_soln;
413, 414 are commented out // since the format to be used for
complex number is unclear
---------------------------------------------------------------------
error:
503,504 and 506 cannot convert ‘std::complex<double>’ to ‘double’
alteration:
456 typedef std::map<int,std::vector<Real> > map_local_soln;
-> typedef std::map<int,std::vector<Number> > map_local_soln;
502 std::vector<Real> vec(3);
-> std::vector<Number> vec(3);
515 - 522 are commented out // since the format to be used for
complex number is unclear
-----------------------------------------------------------------------------
file: src/numerics/numeric_vector.C
error:
221 norm += v(*it) * v(*it); //this is not valid for complex
alteration:
221 norm += std::norm(v(*it));
-----------------------------------------------------------------------------
file: src/numerics/distributed_vector.C
error:
321 this->set(i,fabs(_values[i])); // function 'fabs' is only
defined for double, float, and long double
alterationn:
321 this->set(i,std::abs(_values[i]));
-----------------------------------------------------------------------------
file: src/numerics/laspack_vector.C
error:
262 this->set(i,fabs(_values[i])); // function 'fabs' is only
defined for double, float, and long double
alterationn:
262 this->set(i,std::abs(_values[i]));
-----------------------------------------------------------------------------
also it may be worth mentioning, I always get a warning says that
-----------------------------------------------------------------------------
/libmesh/include/parallel/parallel.h:82: Warning: inline-Function
»Parallel::data_type Parallel::datatype() [with T
= std::complex<double>]« is called, but never defined
-----------------------------------------------------------------------------
Now compilation is no problem any longer, but linking. I got a load of
undefined reference error while building bin/amr. I guess those
functions expect a complex type as parameters, and therefore not
defined. I tried ./configure --disable-amr, --enable-amr=no and etc. to
avoid linking amr, but no luck. Perhaps you can help me out this one,
too :) would be really gradful.
Best regards
Ping
On Fr, 2009-03-06 at 08:31 -0600, Roy Stogner wrote:
> On Fri, 6 Mar 2009, Ping Rong wrote:
>
> > Compiling C++ (in optimized mode) src/mesh/ensight_io.C...
> > src/mesh/ensight_io.C: In member function ‘void
> > EnsightIO::write_scalar_ascii(const std::string&, const
> > std::string&)’:
> > src/mesh/ensight_io.C:407: error: cannot convert
> > ‘std::complex<double>’ to ‘double’ in assignment
>
> Hmm. It's not too uncommon that some jerk adds code to libMesh SVN
> without remembering to check compatibility with --enable-complex, but
> usually I'm the jerk. This one wasn't mine, though, and it's not
> going to be a trivial fix since fixing it properly requires knowing if
> and how the ensight file format handles complex numbers.
>
> > The error in line 503, 504 and 506 are basically the same problem. Could
> > the <Real> type be replaced by <Number> type? Or maybe something else is
> > relying on this.
>
> No, this isn't one of our core classes; you're not going to need it
> unless you explicitly try to read or write a file in that format. In
> fact, if the file's giving you this much trouble I'd say just comment
> the whole thing out (along with ensight_io.h, to avoid linker errors)
> until we get a fix into SVN. Sorry about the trouble.
> ---
> Roy
--
Mit freundlichen Grüßen
M.Sc. Ping Rong
Technische Universität Hamburg-Harburg
Institut für Modellierung und Berechnung
Denickestraße 17
21073 Hamburg
DE 17 Raum 3031
Tel.: ++49 - (0)40 - 42878 2749
Fax: ++49 - (0)40 - 42878 43533

On Fri, 6 Mar 2009, Ping Rong wrote:
> Compiling C++ (in optimized mode) src/mesh/ensight_io.C...
> src/mesh/ensight_io.C: In member function ‘void
> EnsightIO::write_scalar_ascii(const std::string&, const
> std::string&)’:
> src/mesh/ensight_io.C:407: error: cannot convert
> ‘std::complex<double>’ to ‘double’ in assignment
Hmm. It's not too uncommon that some jerk adds code to libMesh SVN
without remembering to check compatibility with --enable-complex, but
usually I'm the jerk. This one wasn't mine, though, and it's not
going to be a trivial fix since fixing it properly requires knowing if
and how the ensight file format handles complex numbers.
> The error in line 503, 504 and 506 are basically the same problem. Could
> the <Real> type be replaced by <Number> type? Or maybe something else is
> relying on this.
No, this isn't one of our core classes; you're not going to need it
unless you explicitly try to read or write a file in that format. In
fact, if the file's giving you this much trouble I'd say just comment
the whole thing out (along with ensight_io.h, to avoid linker errors)
until we get a fix into SVN. Sorry about the trouble.
---
Roy