libmesh-users — General information for libmesh users

[Libmesh-users] Adaptative Mesh Refinement

[Vinicius C. R. <vin...@co...>]

Hi,

I'm trying to perform a h-refinement and a subdomain
evaluation/reevaluation for each refinement iteration, prior to
initializing an EquationSystems object. Currently when call init() I get a
out_of_range exception. I suspect it is thrown from within its node or the
element loop. Everything works fine if only I leave the subdomain
evaluation prior to the refinements (i.e., I comment the subdomain
reevaluation inside the refinement iterations).

The refinement is made by setting the refine flag directly in each selected
element and calling refine_elements() in the MeshRefinement object. I'm not
using any ErrorEstimator because there's no error to be estimated a the
point the refinement is performed.

Can someone point me what might be the wrong assumption(s) I'm making? Is
it possible for a parent element to have childrens in different subdomains?
If I iterate from *MeshBase->active_local_subdomain_elements_begin
([MySubDomain])* through *MeshBase->active_local_subdomain_elements_begin
([MySubDomain])* right after a refinement step will the iterators be
correctly pointing to the new subdomain's range?

I'm using libmesh compiled with Petsc and MPI support, nevertheless I'm
running it in a desktop and with only one node (-np1).

If it helps to contextualize, I'm working on an embedded domain problem.
I'd appreciate any help!

Thanks,
Vinícius da Costa Reis

Master's Candidate
Department of Civil Engineering
COPPE/Federal University of Rio de Janeiro

Re: [Libmesh-users] Variable groups to use blocked matrix (*BAIJ)

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

On Sat, 12 May 2018, Adriano Côrtes wrote:

> If we use a fixed grid with no AMR, then the code works pretty fine,
> but with AMR turned on we are having the same issue as reported by
> David in the mailing list thread pointed by Roy.

> That is, there is a MatSetValue in a position that is not
> preallocated in the BAIJ data structure.  Let us know if we could be
> of any help to fix this issue. The reported error is bellow.

Sorry about the delayed reply.

Is there any possibility you could set up a small test case that
triggers the problem?  I'm running through our usual test suites with
--enable-blocked-storage now to see if I can set it off, but at least
the initial AMR tests are still passing so I'm not sure how likely the
rest are to reproduce the bug.

Thanks,
---
Roy

[Libmesh-users] Solving PDE with subdomains as BCs

[Viviana P. B. <vpa...@uc...>]

Hello,

I am trying to solve the diffusion equation with two types of Dirichlet
BCs. Let's say for example:
df/dt = lapl(f). f(0<x<0.5,y,t) = 1; f(0.5<x<1,y,t) = 0.
I was thinking of taking advantage of subdomains in which subdomain_id = 1
is the bulk. subdomain_id = 2 is surface1. subdomain_id=3 is surface 2.
Once the subdomains are read correctly, I could assemble the RHS
accordingly by using active_local_subdomain_elements_begin( subdomain_id ).
I was able to successfully assemble the corresponding terms using both type
of subdomains but I'm seeing that although it's solved simultaneously, the
solution of one subdomain is not being projected onto the other one. Is
this something I should project manually?

Thank you in advance.

Best,
Viviana.

Re: [Libmesh-users] p-refinement for Lagrange variables?

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

Hi Roy,

Regarding the thread below: I wasn't interested in adaptive p-refinement
per se in this case, instead my motivation was to enable different variable
orders in different regions of a mesh. One use-case for this is to have 2nd
order beam elements and 1st order shell elements in the same mesh, for
example. The beams and shells use the same variables
(u,v,w,theta_x,theta_y,theta_z), but it'd be nice to vary the orders
depending on whether we're in a beam subdomain or a shell subdomain.

It occurred to me that one way to do this more directly than with
p-refinement would be to enable "subdomain based variable order", i.e.
analogously to how we can specify whether a variable is active on a
subdomain, we could also specify its order on a subdomain. Does that seem
like a reasonable idea to you? And what do you think regarding the
implementation, would it be a big undertaking?

Best,
David


On Fri, Mar 9, 2018 at 3:02 PM, David Knezevic <dav...@ak...>
wrote:

> On Fri, Mar 9, 2018 at 2:53 PM, Roy Stogner <roy...@ic...>
> wrote:
>
>>
>> On Fri, 9 Mar 2018, David Knezevic wrote:
>>
>> - If I were to look into adding a special case for non-uniform 1st
>>> to 2nd order refinement for LAGRANGE variables, do you think this
>>> would be of interest to include that in libMesh, or would it be too
>>> specific to include? (I'd like to know if it's potentially of
>>> broader interest before looking further into this.)
>>>
>>
>> Uninteresting to me, but I try like mad to avoid writing
>> LAGRANGE-specific code, so if I want C0 p refinement I just use
>> HIERARCHIC.  Others who've coded themselves too far into a
>> LAGRANGE-only corner might disagree... but from what I've seen the
>> very first way to screw up is to assume that every node has a dof for
>> every variable, and any code like that will *still* be broken if those
>> users have second-order geometric elements (so they can support p=2)
>> but start with p=1.
>>
>> - How complex do you think it would be to add that special case?
>>>
>>
>> Not very.  I personally don't think it's worth it when you'd just end
>> up stuck restricted by p<=2 anyways, but if anyone else disagrees I'd
>> still happily merge their work.  ;-)
>>
>
>
> OK, thanks for that info. I'll think some more about whether this is worth
> working on or not for my use-case.
>
> Thanks,
> David
>

Re: [Libmesh-users] using QGrid rule with very high order

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

On Wed, May 23, 2018 at 12:54 PM, Charles Puelz <cha...@gm...>
wrote:

> thanks for the reply.
>
> this is part of a dynamic quadrature rule, so its not clear what orders I
> may need, but in my current application, they are almost always greater
> than 43.
>

I think it should be possible to use an arbitrary integral value in the
QGrid class.

I opened Github issue https://github.com/libMesh/libmesh/issues/1711 to
keep the discussion going.

-- 
John

Re: [Libmesh-users] Use of comm.barrier()

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

On Wed, May 23, 2018 at 6:37 PM, Viviana Palacio Betancur <
vpa...@uc...> wrote:

> Hello,
>
> I have a correct solution of my problem and now I want to run in parallel.
> I'm using all the correct iterators and the program compiles and runs,
> except now the solution is completely different and inconsistent to the
> problem run in serial. I'm using the SuperLU solver.
>
> I suspect that I'm making some mistakes when running in parallel, since it
> might be the case that the information is not being shared simultaneously?
> Is it necessary to impose communication barriers and if so in which steps
> of the solution should this be implemented?
>

In general it is not necessary to impose communication barriers in order to
get library code working in parallel, I can't speak for any parallel code
you may have written yourself.

The usual issue is that you are using an active_element_iterator for
assembly, which works fine in serial, but obviously breaks when you try to
switch to parallel because you should be using
active_local_element_iterators.

Another possibility is that you are missing a system.update() call
somewhere. This can cause issues if you are using current_local_solution
values before they have been updated with values from the parallel solution
vector, but isn't very common.

Without seeing your code, it's hard to say what the issue might be...

-- 
John

[Libmesh-users] Matrix Preallocation for continuous FE type (LAGRANGE) and discontinuous FE type (L2_LAGRANGE)

[<li...@si...>]

Dear developers,


I want to couple two variables into one system, one with the continuous FE type and the other with the discontinuous FE type.


I tried 
            system.add_variable ("p", FIRST, LAGRANGE);
            system.add_variable ("s", FIRST, L2_LAGRANGE);


When inserting the neighbor matrix "Ken" for the second variable into the global stiffness matrix as in "example/miscellaneous_ex5", PETSc error reports that 
there is no MatrixPreallocation for inserting "Ken". 


But if I used L2_LAGRANGE for both variables, 
            system.add_variable ("p", FIRST, L2_LAGRANGE);
            system.add_variable ("s", FIRST, L2_LAGRANGE);
There is no error report. I believe the sparsity_patterns are different between LAGRANGE and L2_LAGRANGE. 


But right now I need to use different FE types for the two variable respectively.
I don't know why it outputs error when I do the inserting of "Ken" ONLY for the second variable, which is already set to be L2_LAGRANGE. 
It seems that as long as one variable in the system is set to LAGRANGE, the whole Dof_map will build a sparsity_pattern based on the LAGRANGE type, even the other variable is set to L2_LAGRANGE type. It is hard for me to find such setting in the libmesh source codes such as dof_map.C.




Best,
Li Luo


















-----Original Messages-----
From:"John Peterson" <jwp...@gm...>
Sent Time:2018-05-17 00:19:38 (Thursday)
To: "li.luo" <li...@si...>
Cc: libmesh-users <lib...@li...>
Subject: Re: [Libmesh-users] The use of MONOMIAL/XYZ FE_type






On Wed, May 16, 2018 at 6:53 AM, <li...@si...> wrote:
Dear developers,




I tried using the MONOMIAL/XYZ FE_type as in the miscellaneous_ex5 example.


After simplification under the choice of some parameters, my problem reduces to a Laplace equation with Dirichlet boundary condition.
If a default FE_type is used (system.add_variable ("p", FIRST);), the Laplace problem outputs a right result,
but if a MONOMIAL/XYZ FE_type is used instead (i.e. system.add_variable ("p", FIRST, XYZ)), 
the interior part is wrong, it isn't changed from its initial solution value (I set zero for the initial value),
only the boundary is right due to the impose of Dirichlet condition.
I have implemented the contribution of each element boundary flux as that in the miscellaneous_ex5 example.
I don't know why it fails in the interior part.



I'm not sure why XYZ would not work. It's possible that the solution is actually OK but the error arises in the ExodusII::write_discontinuous_exodusII() output? 


I don't think we have much testing of the XYZ FEType in libmesh or in apps which use libmesh, so it's possible that i regressed.




 
2.  Another question is, if I want to couple two variables into one system, one with a default FE_type and the other with the MONOMIAL/XYZ FE_type, is that possible?
I tried and found that when inserting the contribution of neighbor values into the matrix (as "Ken, Kne" in miscellaneous_ex5),
PETSc error reports that
---------------------------------------------------------------
Argument out of range!
New nonzero at (13,33164) caused a malloc!
MatSetValues_MPIAIJ() line 572 in src/mat/impls/aij/mpi/mpiaij.c
MatSetValues() line 1106 in src/mat/interface/matrix.c
add_matrix() line 699 in "unknowndirectory/"petsc_matrix.C
---------------------------------------------------------------


It seems that different FE_types lead to different matrix preallocations,
the MONOMIAL/XYZ FE_type allows the insert of Ken, Kne but the default FE_type does not.
Can I use the same matrix preallocation for the two variables to allow such inserting?



And you are using libmesh master? This is interesting, as we are discussing a similar error right now over on https://github.com/libMesh/libmesh/issues/1697 that might be related, although I suspect yours might have more to do with DG sparsity patterns.


--

John

--
Shenzhen Institutes of Advanced Technology, Chinese Academy of Sciences 1068 Xueyuan Avenue, Shenzhen University Town, Shenzhen, P. R. China Tel: +86 13691688103 Mail: li...@si... https://sites.google.com/site/rolyliluo/

[Libmesh-users] Use of comm.barrier()

[Viviana P. B. <vpa...@uc...>]

Hello,

I have a correct solution of my problem and now I want to run in parallel.
I'm using all the correct iterators and the program compiles and runs,
except now the solution is completely different and inconsistent to the
problem run in serial. I'm using the SuperLU solver.

I suspect that I'm making some mistakes when running in parallel, since it
might be the case that the information is not being shared simultaneously?
Is it necessary to impose communication barriers and if so in which steps
of the solution should this be implemented?

Thank you. Best,

Viviana.

Re: [Libmesh-users] using QGrid rule with very high order

[Charles P. <cha...@gm...>]

thanks for the reply.

this is part of a dynamic quadrature rule, so its not clear what orders I
may need, but in my current application, they are almost always greater
than 43.

On Wed, May 23, 2018 at 2:51 PM, John Peterson <jwp...@gm...> wrote:

>
>
> On Wed, May 23, 2018 at 11:10 AM, Charles Puelz <cha...@gm...>
> wrote:
>
>> Hi everyone,
>>
>> I'd like to use the QGrid object as a quadrature rule, and use an order
>> that's potentially
>> much bigger than what is defined in the Order enum.  For example, I'd like
>> to do:
>>
>> Order foo = static_cast<Order>(100)
>>
>> but I am concerned about undefined behavior with this cast.  Is there a
>> good way
>> to use QGrid with an order greater than 43?
>>
>
> The easiest thing would probably be to just add your own enumerations up
> to whatever order you need. The list can be sparse, i.e. you don't have to
> list everything between 43 and 100.
>
> --
> John
>



-- 
Charles Puelz
Postdoc
UNC Math Department
cpuelz.github.io

Re: [Libmesh-users] using QGrid rule with very high order

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

On Wed, May 23, 2018 at 11:10 AM, Charles Puelz <cha...@gm...>
wrote:

> Hi everyone,
>
> I'd like to use the QGrid object as a quadrature rule, and use an order
> that's potentially
> much bigger than what is defined in the Order enum.  For example, I'd like
> to do:
>
> Order foo = static_cast<Order>(100)
>
> but I am concerned about undefined behavior with this cast.  Is there a
> good way
> to use QGrid with an order greater than 43?
>

The easiest thing would probably be to just add your own enumerations up to
whatever order you need. The list can be sparse, i.e. you don't have to
list everything between 43 and 100.

-- 
John

[Libmesh-users] using QGrid rule with very high order

[Charles P. <cha...@gm...>]

Hi everyone,

I'd like to use the QGrid object as a quadrature rule, and use an order
that's potentially
much bigger than what is defined in the Order enum.  For example, I'd like
to do:

Order foo = static_cast<Order>(100)

but I am concerned about undefined behavior with this cast.  Is there a
good way
to use QGrid with an order greater than 43?

Thanks!

-- Charles

-- 
Charles Puelz
Postdoc
UNC Math Department
cpuelz.github.io

[Libmesh-users] flat subdomain boundary

[Chang, H. <mo...@nf...>]

Hi,



Now I try to solve Poisson's equation for the electrostatic field and continuity equations for plasma fluids in the same mesh.



In this case, dielectric region has to be included in Poisson solve, but not in continuity equations.



In LIBMESH example, I found that subdomain method make it possible to determine whether included or not.



I'd like to ask about how can I make flat the uneven boundaries between two different domains.



Not only is it not good to look at but it also presents difficulties in correcting errors in boundary conditions.

Re: [Libmesh-users] Connection between a 2D and 3D mesh

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

On Thu, 17 May 2018, Viviana Palacio Betancur wrote:

> Thank you for your response! I was able to find the parent element. Is
> there a way I could obtain the nodal information as well?

Elem::node_ref, node_ptr, or node_id will get you a reference to,
pointer to, or id() of a node indexed by local node number.

If you're looking specifically for the mapping between local node
numbers on a surface element vs on the interior_parent() then there's
no single API for that, but you can call node_id or node_ptr on one
and then look up the result with local_node or get_node_index on the
other.
---
Roy

Re: [Libmesh-users] Connection between a 2D and 3D mesh

[Viviana P. B. <vpa...@uc...>]

Thank you for your response! I was able to find the parent element. Is
there a way I could obtain the nodal information as well?

Viviana Palacio Betancur
PhD Student | 2016 Cohort
de Pablo Group
Institute for Molecular Engineering
University of Chicago

On Thu, May 17, 2018 at 4:34 PM, John Peterson <jwp...@gm...> wrote:

>
>
> On Thu, May 17, 2018 at 2:58 PM, Viviana Palacio Betancur <
> vpa...@uc...> wrote:
>
>> Hello,
>>
>> Currently I'm working on calculating surface normals for a 3D mesh. The
>> way
>> I have set up this is by extracting the surface mesh by using
>> mesh->get_boundary_info().sync(mesh_surf);
>> This gives me a 2D mesh that I can use for calculating element normals.
>> However, I haven't been able to establish the connection between the
>> element ids of the surface mesh to the ids of the original mesh.
>>
>
> The elements in mesh_surf will have an interior_parent() which points to
> the higher-dimensional element they came from.
>
>
> --
> John
>

Re: [Libmesh-users] Connection between a 2D and 3D mesh

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

On Thu, May 17, 2018 at 2:58 PM, Viviana Palacio Betancur <
vpa...@uc...> wrote:

> Hello,
>
> Currently I'm working on calculating surface normals for a 3D mesh. The way
> I have set up this is by extracting the surface mesh by using
> mesh->get_boundary_info().sync(mesh_surf);
> This gives me a 2D mesh that I can use for calculating element normals.
> However, I haven't been able to establish the connection between the
> element ids of the surface mesh to the ids of the original mesh.
>

The elements in mesh_surf will have an interior_parent() which points to
the higher-dimensional element they came from.


-- 
John

[Libmesh-users] Connection between a 2D and 3D mesh

[Viviana P. B. <vpa...@uc...>]

Hello,

Currently I'm working on calculating surface normals for a 3D mesh. The way
I have set up this is by extracting the surface mesh by using
mesh->get_boundary_info().sync(mesh_surf);
This gives me a 2D mesh that I can use for calculating element normals.
However, I haven't been able to establish the connection between the
element ids of the surface mesh to the ids of the original mesh.

Is there an existing id list from when I first performed
get_boundary_info().

Thank you. Best,
Viviana.


Viviana Palacio Betancur

Re: [Libmesh-users] The use of MONOMIAL/XYZ FE_type

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

On Thu, May 17, 2018 at 2:13 AM, <li...@si...> wrote:

>
> Thank you for the reply.
> However, I still don't know how to fix it. Any hints?
>

Nothing beyond good old fashioned debugging. As I suggested, you could
start by checking whether the problem lies in the FE formulation/solution
itself or is an error within the Exodus writer.


> I am not using libmesh master.
>

If you plan to do serious debugging, I highly recommend working with
libmesh master. Otherwise, anything you fix may not be applicable/mergeable
to the current version of the code.

-- 
John

Re: [Libmesh-users] The use of MONOMIAL/XYZ FE_type

[<li...@si...>]

Thank you for the reply.
However, I still don't know how to fix it. Any hints?

I am not using libmesh master.


Regards,
Li Luo



-----Original Messages-----
From:"John Peterson" <jwp...@gm...>
Sent Time:2018-05-17 00:19:38 (Thursday)
To: "li.luo" <li...@si...>
Cc: libmesh-users <lib...@li...>
Subject: Re: [Libmesh-users] The use of MONOMIAL/XYZ FE_type






On Wed, May 16, 2018 at 6:53 AM, <li...@si...> wrote:
Dear developers,




I tried using the MONOMIAL/XYZ FE_type as in the miscellaneous_ex5 example.


After simplification under the choice of some parameters, my problem reduces to a Laplace equation with Dirichlet boundary condition.
If a default FE_type is used (system.add_variable ("p", FIRST);), the Laplace problem outputs a right result,
but if a MONOMIAL/XYZ FE_type is used instead (i.e. system.add_variable ("p", FIRST, XYZ)), 
the interior part is wrong, it isn't changed from its initial solution value (I set zero for the initial value),
only the boundary is right due to the impose of Dirichlet condition.
I have implemented the contribution of each element boundary flux as that in the miscellaneous_ex5 example.
I don't know why it fails in the interior part.



I'm not sure why XYZ would not work. It's possible that the solution is actually OK but the error arises in the ExodusII::write_discontinuous_exodusII() output? 


I don't think we have much testing of the XYZ FEType in libmesh or in apps which use libmesh, so it's possible that i regressed.




 
2.  Another question is, if I want to couple two variables into one system, one with a default FE_type and the other with the MONOMIAL/XYZ FE_type, is that possible?
I tried and found that when inserting the contribution of neighbor values into the matrix (as "Ken, Kne" in miscellaneous_ex5),
PETSc error reports that
---------------------------------------------------------------
Argument out of range!
New nonzero at (13,33164) caused a malloc!
MatSetValues_MPIAIJ() line 572 in src/mat/impls/aij/mpi/mpiaij.c
MatSetValues() line 1106 in src/mat/interface/matrix.c
add_matrix() line 699 in "unknowndirectory/"petsc_matrix.C
---------------------------------------------------------------


It seems that different FE_types lead to different matrix preallocations,
the MONOMIAL/XYZ FE_type allows the insert of Ken, Kne but the default FE_type does not.
Can I use the same matrix preallocation for the two variables to allow such inserting?



And you are using libmesh master? This is interesting, as we are discussing a similar error right now over on https://github.com/libMesh/libmesh/issues/1697 that might be related, although I suspect yours might have more to do with DG sparsity patterns.


--

John

--
Shenzhen Institutes of Advanced Technology, Chinese Academy of Sciences 1068 Xueyuan Avenue, Shenzhen University Town, Shenzhen, P. R. China Tel: +86 13691688103 Mail: li...@si... https://sites.google.com/site/rolyliluo/

Re: [Libmesh-users] The use of MONOMIAL/XYZ FE_type

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

On Wed, May 16, 2018 at 6:53 AM, <li...@si...> wrote:

> Dear developers,
>
>
>
>
> I tried using the MONOMIAL/XYZ FE_type as in the miscellaneous_ex5 example.
>
>
> After simplification under the choice of some parameters, my problem
> reduces to a Laplace equation with Dirichlet boundary condition.
> If a default FE_type is used (system.add_variable ("p", FIRST);), the
> Laplace problem outputs a right result,
> but if a MONOMIAL/XYZ FE_type is used instead (i.e. system.add_variable
> ("p", FIRST, XYZ)),
> the interior part is wrong, it isn't changed from its initial solution
> value (I set zero for the initial value),
> only the boundary is right due to the impose of Dirichlet condition.
> I have implemented the contribution of each element boundary flux as that
> in the miscellaneous_ex5 example.
> I don't know why it fails in the interior part.
>

I'm not sure why XYZ would not work. It's possible that the solution is
actually OK but the error arises in the
ExodusII::write_discontinuous_exodusII() output?

I don't think we have much testing of the XYZ FEType in libmesh or in apps
which use libmesh, so it's possible that i regressed.




> 2.  Another question is, if I want to couple two variables into one
> system, one with a default FE_type and the other with the MONOMIAL/XYZ
> FE_type, is that possible?
> I tried and found that when inserting the contribution of neighbor values
> into the matrix (as "Ken, Kne" in miscellaneous_ex5),
> PETSc error reports that
> ---------------------------------------------------------------
> Argument out of range!
> New nonzero at (13,33164) caused a malloc!
> MatSetValues_MPIAIJ() line 572 in src/mat/impls/aij/mpi/mpiaij.c
> MatSetValues() line 1106 in src/mat/interface/matrix.c
> add_matrix() line 699 in "unknowndirectory/"petsc_matrix.C
> ---------------------------------------------------------------
>
>
> It seems that different FE_types lead to different matrix preallocations,
> the MONOMIAL/XYZ FE_type allows the insert of Ken, Kne but the default
> FE_type does not.
> Can I use the same matrix preallocation for the two variables to allow
> such inserting?
>

And you are using libmesh master? This is interesting, as we are discussing
a similar error right now over on
https://github.com/libMesh/libmesh/issues/1697 that might be related,
although I suspect yours might have more to do with DG sparsity patterns.

-- 
John

[Libmesh-users] The use of MONOMIAL/XYZ FE_type

[<li...@si...>]

Dear developers,




I tried using the MONOMIAL/XYZ FE_type as in the miscellaneous_ex5 example.


After simplification under the choice of some parameters, my problem reduces to a Laplace equation with Dirichlet boundary condition.
If a default FE_type is used (system.add_variable ("p", FIRST);), the Laplace problem outputs a right result, 
but if a MONOMIAL/XYZ FE_type is used instead (i.e. system.add_variable ("p", FIRST, XYZ)),  
the interior part is wrong, it isn't changed from its initial solution value (I set zero for the initial value), 
only the boundary is right due to the impose of Dirichlet condition. 
I have implemented the contribution of each element boundary flux as that in the miscellaneous_ex5 example. 
I don't know why it fails in the interior part.




2.  Another question is, if I want to couple two variables into one system, one with a default FE_type and the other with the MONOMIAL/XYZ FE_type, is that possible?
I tried and found that when inserting the contribution of neighbor values into the matrix (as "Ken, Kne" in miscellaneous_ex5), 
PETSc error reports that 
---------------------------------------------------------------
Argument out of range!
New nonzero at (13,33164) caused a malloc!
MatSetValues_MPIAIJ() line 572 in src/mat/impls/aij/mpi/mpiaij.c
MatSetValues() line 1106 in src/mat/interface/matrix.c
add_matrix() line 699 in "unknowndirectory/"petsc_matrix.C
---------------------------------------------------------------


It seems that different FE_types lead to different matrix preallocations,
the MONOMIAL/XYZ FE_type allows the insert of Ken, Kne but the default FE_type does not.
Can I use the same matrix preallocation for the two variables to allow such inserting?


Thank you!


Regards,
Li Luo

Re: [Libmesh-users] normal vectors on surface meshes

[Griffith, B. E. <bo...@em...>]

On May 15, 2018, at 10:23 AM, John Peterson <jwp...@gm...<mailto:jwp...@gm...>> wrote:



On Tue, May 15, 2018 at 8:17 AM, Griffith, Boyce Eugene <bo...@em...<mailto:bo...@em...>> wrote:


On May 15, 2018, at 10:13 AM, John Peterson <jwp...@gm...<mailto:jwp...@gm...>> wrote:



On Tue, May 15, 2018 at 5:37 AM, Griffith, Boyce Eugene <bo...@em...<mailto:bo...@em...>> wrote:
Folks --

Is there an easy way in the library to get normal vectors on surface meshes (e.g. 2D elements in 3D space, 1D elements in 2D space, etc)?  It seems like most of the support for this assumes that these are only needed for doing surface integrals on volumetric meshes, but I am sure I am overlooking something.

The "simplest" way is to build an FE object and only pre-request the normals by calling:

const std::vector<Point> & face_normals = fe_face->get_normals();

I may be making a blunder here, but this doesn't seem to work if the dimension of the quadrature rule is equal to the dimension of the mesh --- the normals do not appear to be computed.


Hmm... yes generally you would use a dim-1 dimensional quadrature rule in this context, e.g. the following (from adaptivity_ex4) should work:

std::unique_ptr<FEBase> fe_face (FEBase::build(dim, fe_type));
std::unique_ptr<QBase> qface(fe_type.default_quadrature_rule(dim-1));
fe_face->attach_quadrature_rule (qface.get());
const std::vector<Point> & face_normals = fe_face->get_normals();
fe_face->reinit(elem, s);

I need to integrate on a surface mesh (for instance, a boundary mesh extracted from a volumetric mesh --- although it could simply be a surface that is not derived from a volumetric mesh) and to get the normals --- something like:

std::unique_ptr<FEBase> fe_surface(FEBase::build(dim, fe_type));
std::unique_ptr<QBase> quad_surface(fe_type.default_quadrature_rule(dim));
fe_surface->attach_quadrature_rule (quad_surface.get());
const std::vector<Point> & face_normals = fe_surface->get_normals();
fe_surface->reinit(elem);

Again, I do know that I can do this myself --- I was just going to see if the library already does it.

where "s" is the numerical index of a side.

--
John

Re: [Libmesh-users] normal vectors on surface meshes

[Paul T. B. <ptb...@gm...>]

Not to butt-in, but just to confirm, I also compute this "manually" as John
suggests for things like normal pressure on membranes, etc. I don't believe
there's a canned function for it.

On Tue, May 15, 2018 at 10:38 AM, John Peterson <jwp...@gm...>
wrote:

> On Tue, May 15, 2018 at 8:31 AM, Griffith, Boyce Eugene <
> bo...@em...> wrote:
>
> >
> >
> > On May 15, 2018, at 10:23 AM, John Peterson <jwp...@gm...>
> wrote:
> >
> >
> >
> > On Tue, May 15, 2018 at 8:17 AM, Griffith, Boyce Eugene <
> > bo...@em...> wrote:
> >
> >>
> >>
> >> On May 15, 2018, at 10:13 AM, John Peterson <jwp...@gm...>
> wrote:
> >>
> >>
> >>
> >> On Tue, May 15, 2018 at 5:37 AM, Griffith, Boyce Eugene <
> >> bo...@em...> wrote:
> >>
> >>> Folks --
> >>>
> >>> Is there an easy way in the library to get normal vectors on surface
> >>> meshes (e.g. 2D elements in 3D space, 1D elements in 2D space, etc)?
> It
> >>> seems like most of the support for this assumes that these are only
> needed
> >>> for doing surface integrals on volumetric meshes, but I am sure I am
> >>> overlooking something.
> >>>
> >>
> >> The "simplest" way is to build an FE object and only pre-request the
> >> normals by calling:
> >>
> >> const std::vector<Point> & face_normals = fe_face->get_normals();
> >>
> >>
> >> I may be making a blunder here, but this doesn't seem to work if the
> >> dimension of the quadrature rule is equal to the dimension of the mesh
> ---
> >> the normals do not appear to be computed.
> >>
> >
> >
> > Hmm... yes generally you would use a dim-1 dimensional quadrature rule in
> > this context, e.g. the following (from adaptivity_ex4) should work:
> >
> > std::unique_ptr<FEBase> fe_face (FEBase::build(dim, fe_type));
> > std::unique_ptr<QBase> qface(fe_type.default_quadrature_rule(dim-1));
> > fe_face->attach_quadrature_rule (qface.get());
> > const std::vector<Point> & face_normals = fe_face->get_normals();
> > fe_face->reinit(elem, s);
> >
> >
> > I need to integrate on a surface mesh (for instance, a boundary mesh
> > extracted from a volumetric mesh --- although it could simply be a
> surface
> > that is not derived from a volumetric mesh) and to get the normals ---
> > something like:
> >
> > std::unique_ptr<FEBase> fe_surface(FEBase::build(dim, fe_type));
> > std::unique_ptr<QBase> quad_surface(fe_type.default_
> quadrature_rule(dim));
> > fe_surface->attach_quadrature_rule (quad_surface.get());
> > const std::vector<Point> & face_normals = fe_surface->get_normals();
> > fe_surface->reinit(elem);
> >
>
>
> I see now what you are saying. Yeah, the normals are only computed during
> *side* reinits, since they are outward normals relative to
> higher-dimensional element.
>
> So, for a 2D element I would compute them manually by computing dx/d(xi)
> \cross dx/d(eta), (those two quantities are available from the FE object to
> help you).
>
> --
> John
> ------------------------------------------------------------
> ------------------
> Check out the vibrant tech community on one of the world's most
> engaging tech sites, Slashdot.org! http://sdm.link/slashdot
> _______________________________________________
> Libmesh-users mailing list
> Lib...@li...
> https://lists.sourceforge.net/lists/listinfo/libmesh-users
>

Re: [Libmesh-users] normal vectors on surface meshes

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

On Tue, May 15, 2018 at 8:31 AM, Griffith, Boyce Eugene <
bo...@em...> wrote:

>
>
> On May 15, 2018, at 10:23 AM, John Peterson <jwp...@gm...> wrote:
>
>
>
> On Tue, May 15, 2018 at 8:17 AM, Griffith, Boyce Eugene <
> bo...@em...> wrote:
>
>>
>>
>> On May 15, 2018, at 10:13 AM, John Peterson <jwp...@gm...> wrote:
>>
>>
>>
>> On Tue, May 15, 2018 at 5:37 AM, Griffith, Boyce Eugene <
>> bo...@em...> wrote:
>>
>>> Folks --
>>>
>>> Is there an easy way in the library to get normal vectors on surface
>>> meshes (e.g. 2D elements in 3D space, 1D elements in 2D space, etc)?  It
>>> seems like most of the support for this assumes that these are only needed
>>> for doing surface integrals on volumetric meshes, but I am sure I am
>>> overlooking something.
>>>
>>
>> The "simplest" way is to build an FE object and only pre-request the
>> normals by calling:
>>
>> const std::vector<Point> & face_normals = fe_face->get_normals();
>>
>>
>> I may be making a blunder here, but this doesn't seem to work if the
>> dimension of the quadrature rule is equal to the dimension of the mesh ---
>> the normals do not appear to be computed.
>>
>
>
> Hmm... yes generally you would use a dim-1 dimensional quadrature rule in
> this context, e.g. the following (from adaptivity_ex4) should work:
>
> std::unique_ptr<FEBase> fe_face (FEBase::build(dim, fe_type));
> std::unique_ptr<QBase> qface(fe_type.default_quadrature_rule(dim-1));
> fe_face->attach_quadrature_rule (qface.get());
> const std::vector<Point> & face_normals = fe_face->get_normals();
> fe_face->reinit(elem, s);
>
>
> I need to integrate on a surface mesh (for instance, a boundary mesh
> extracted from a volumetric mesh --- although it could simply be a surface
> that is not derived from a volumetric mesh) and to get the normals ---
> something like:
>
> std::unique_ptr<FEBase> fe_surface(FEBase::build(dim, fe_type));
> std::unique_ptr<QBase> quad_surface(fe_type.default_quadrature_rule(dim));
> fe_surface->attach_quadrature_rule (quad_surface.get());
> const std::vector<Point> & face_normals = fe_surface->get_normals();
> fe_surface->reinit(elem);
>


I see now what you are saying. Yeah, the normals are only computed during
*side* reinits, since they are outward normals relative to
higher-dimensional element.

So, for a 2D element I would compute them manually by computing dx/d(xi)
\cross dx/d(eta), (those two quantities are available from the FE object to
help you).

-- 
John

Re: [Libmesh-users] normal vectors on surface meshes

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

On Tue, May 15, 2018 at 8:17 AM, Griffith, Boyce Eugene <
bo...@em...> wrote:

>
>
> On May 15, 2018, at 10:13 AM, John Peterson <jwp...@gm...> wrote:
>
>
>
> On Tue, May 15, 2018 at 5:37 AM, Griffith, Boyce Eugene <
> bo...@em...> wrote:
>
>> Folks --
>>
>> Is there an easy way in the library to get normal vectors on surface
>> meshes (e.g. 2D elements in 3D space, 1D elements in 2D space, etc)?  It
>> seems like most of the support for this assumes that these are only needed
>> for doing surface integrals on volumetric meshes, but I am sure I am
>> overlooking something.
>>
>
> The "simplest" way is to build an FE object and only pre-request the
> normals by calling:
>
> const std::vector<Point> & face_normals = fe_face->get_normals();
>
>
> I may be making a blunder here, but this doesn't seem to work if the
> dimension of the quadrature rule is equal to the dimension of the mesh ---
> the normals do not appear to be computed.
>


Hmm... yes generally you would use a dim-1 dimensional quadrature rule in
this context, e.g. the following (from adaptivity_ex4) should work:

std::unique_ptr<FEBase> fe_face (FEBase::build(dim, fe_type));
std::unique_ptr<QBase> qface(fe_type.default_quadrature_rule(dim-1));
fe_face->attach_quadrature_rule (qface.get());
const std::vector<Point> & face_normals = fe_face->get_normals();
fe_face->reinit(elem, s);

where "s" is the numerical index of a side.

-- 
John

Re: [Libmesh-users] normal vectors on surface meshes

[Griffith, B. E. <bo...@em...>]

On May 15, 2018, at 10:13 AM, John Peterson <jwp...@gm...<mailto:jwp...@gm...>> wrote:



On Tue, May 15, 2018 at 5:37 AM, Griffith, Boyce Eugene <bo...@em...<mailto:bo...@em...>> wrote:
Folks --

Is there an easy way in the library to get normal vectors on surface meshes (e.g. 2D elements in 3D space, 1D elements in 2D space, etc)?  It seems like most of the support for this assumes that these are only needed for doing surface integrals on volumetric meshes, but I am sure I am overlooking something.

The "simplest" way is to build an FE object and only pre-request the normals by calling:

const std::vector<Point> & face_normals = fe_face->get_normals();

I may be making a blunder here, but this doesn't seem to work if the dimension of the quadrature rule is equal to the dimension of the mesh --- the normals do not appear to be computed.

(an example can be seen in adaptivity_ex4). The normals could of course be computed using only geometric information, but this is not currently implemented.

--
John