libmesh-users — General information for libmesh users

[Libmesh-users] [RB] Assembly for a Interior (neighbor) side between elements

[Shinseong K. <ss...@pu...>]

<mailplughead>Hello, all.I triedto definean output in a simply supported 3-D beam problem.My output is an average displacement in the middle of the beam.For simplicity, I only used two brick elements, i.e, ■|■, where"■" is an element, and"|" is a location I want to compute the average output.To do this, I utilized the codes in RB example 5.First, Ideclared "ThetaOutput" (parameter-dependent term) and alsodeclared "AssemblyOutput" (parameter-independent term) as follows:=======================================voidAssemblyOutput::boundary_assembly(FEMContext & c)}=======================================As you can see, The above code is very similar to"Assemblyp1" in "assembly.C" ofRB example 5.However, I could not obtain the correct output because the above code cannot deal with interior side between elements.In other words, theboundary information of FEMContext only include outer side.I tried to find another way, but I could not find the right solution.Therefore,I'd like your help.To summarize my questions areas follows:1. Is there any way to assemble an interior (neighbor) sidebetween elements?2. How to compute anaverage displacement for the middleof a beam in terms of the RB codes?Thank you.Best regards,Kang------------------------------------------------------------
ShinseongKang
GraduateStudent
PusanNationalUniversity, South KoreaTel.: +82-051-510-3052
H.P.:+82-010-9770-6595
E-mail:ss...@pu...

------------------------------------------------------------</mailplughtml>

Re: [Libmesh-users] Seg Fault in mesh clone

[Stogner, R. H <roy...@ic...>]

On Sun, 6 Oct 2019, Renato Poli wrote:

> I am trying to clone a mesh so I can remesh the clone and transfer the
> solution from the original mesh to the new one. (is that the best way to
> do??)

It depends on how extensive the remeshing you plan is.  If you can
create the new mesh via isotropic Adaptive Mesh Refinement /
Coarsening of the old one, then just using our built-in AMR/C
facilities is the best way to go.  If you need to do anisotropic
refinement but you're doing something like edge splitting, writing
your own interpolation code to go along with that might be easy
enough and would definitely be most efficient.  But if you're creating
a drastically different mesh then yeah, a clone and transfer is
probably best.

> I am seeing a segmentation fault while calling 'ghosting functors' inside
> 'prepare_for_use'.
> (I do have 'coupling functors' attached to the original System.)
> I believe it should work ... Do you agree?

I agree that it should work, in the normative if not the positive
sense of "should"... but I can see why it might not work the way you
expect; I might not have been getting clone() right with
GhostingFunctor objects, since we don't have much test coverage using
both.

So let's see... when we clone a ReplicatedMesh or a DistributedMesh, a
big hunk of that work is getting done by the (defaulted)
UnstructuredMesh copy constructor, which will in turn call the
MeshBase copy constructor.  The _default_ghosting gets created anew,
but the _ghosting_functors set just gets a shallow copy.


> 302           gf->mesh_reinit();

This is presumably a call to a dangling pointer, then?  You add a
ghosting functor to the old mesh, but you (or the old mesh's
_shared_functors map) have it get deleted afterward, and then we get a
segfault when the new mesh tries to still call it?  Could that be
what's going on?
---
Roy

[Libmesh-users] Seg Fault in mesh clone

[Renato P. <re...@gm...>]

Hi All,

I am trying to clone a mesh so I can remesh the clone and transfer the
solution from the original mesh to the new one. (is that the best way to
do??)
I am seeing a segmentation fault while calling 'ghosting functors' inside
'prepare_for_use'.
(I do have 'coupling functors' attached to the original System.)
I believe it should work ... Do you agree?

Please take a look at the error below.

>> std::unique_ptr<MeshBase> mesh = _old_es->get_mesh().clone();

------
Program received signal SIGSEGV, Segmentation fault.
0x00007ffff2adbabc in libMesh::MeshBase::prepare_for_use
(this=0x555555c77290, skip_renumber_nodes_and_elements=false,
skip_find_neighbors=true) at src/mesh/mesh_base.C:302
302           gf->mesh_reinit();
(gdb) bt
#0  0x00007ffff2adbabc in libMesh::MeshBase::prepare_for_use
(this=0x555555c77290, skip_renumber_nodes_and_elements=false,
skip_find_neighbors=true) at src/mesh/mesh_base.C:302
#1  0x00007ffff2d6b075 in
libMesh::UnstructuredMesh::copy_nodes_and_elements (this=0x555555c77290,
other_mesh=..., skip_find_neighbors=true, element_id_offset=0,
node_id_offset=0)
    at src/mesh/unstructured_mesh.C:218
#2  0x00007ffff2d46825 in libMesh::ReplicatedMesh::ReplicatedMesh
(this=0x555555c77290, other_mesh=...) at src/mesh/replicated_mesh.C:110
#3  0x0000555555667104 in std::make_unique<libMesh::ReplicatedMesh,
libMesh::ReplicatedMesh const&> () at
/usr/include/c++/7/bits/unique_ptr.h:825
#4  libMesh::ReplicatedMesh::clone (this=0x7fffffffcb20) at
/usr/local/include/libmesh/replicated_mesh.h:86
#5  chimas_clib::Chimas::move_tip (this=0x7fffffffd000, dist=5) at
../../../src/clib/chimas_clib/Chimas.cxx:427
-----

rgds,
Renato

Re: [Libmesh-users] Adapting systems_of_equations_ex7.C to heterogeneous Dirichlet BC

[Yuxiang W. <yw...@vi...>]

Hi Roy,

Fantastic! Thank you so much for your help. It works great for us now.

Best,
Shawn

On Tue, Oct 1, 2019 at 1:41 PM Stogner, Roy H <roy...@ic...>
wrote:

>
> On Mon, 30 Sep 2019, Yuxiang Wang wrote:
>
> > I am trying to adapt the systems_of_equations_ex7.C to heterogeneous
> > Dirichlet BC, and I seem to hit a few challenges. Could you please help
> me
> > out here?
> >
> > 1) My first attempt is simply adding the heterogeneous BC like a linear
> > problem and used
> > `dof_map.heterogenously_constrain_element_matrix_and_vector (Ke, Fe,
> > dof_indices)` instead of the `dof_map.constrain_element_matrix_and_vector
> > (Ke, Fe, dof_indices)` in the jacobian() function. Note that I did NOT
> > have `system.rhs->add_vector (Fe, dof_indices)` in the jacobian()
> function.
> > The solution seems to be slightly wrong - in the result file, some nodes
> > did not have the right prescribed displacements, although their
> neighboring
> > nodes seemed correct. Then I added the `system.rhs->add_vector (Fe,
> > dof_indices)` in the jacobian() function then the results diverged.
>
> Treating the heterogeneous BC like a linear problem is tricky if you
> don't have control over what the solver's doing under the hood.
>
> What we tend to do now in that case is a cheat: treat the constrained
> DoFs as "mutable" in the C++ sense, do enforce_constraints_exactly()
> before each residual or Jacobian evaluation, do homogeneous assembly
> during the evaluations, then do enforce_constraints_exactly() one
> final time after the solver returns.  See petsc_diff_solver.C for
> example.
>
> > 2) I noticed that the BCs have to be applied before the es.init().
>
> > However, for NonlinearImplicitSystem we may want to divide the
> > solution into multiple steps and we'd then need to update the
> > heterogeneous BC per step.  Would this be possible?
>
> If you have nonlinear (or time-varying) heterogeneous Dirichlet BCs,
> System::reinit_constraints() can be used to recalculate the constraint
> equations from them when you need.
> ---
> Roy
>


-- 
Yuxiang "Shawn" Wang, PhD
yw...@vi...
+1 (434) 284-0836

Re: [Libmesh-users] Exchanging information between meshes

[Renato P. <re...@gm...>]

Thanks!

On Tue, Oct 1, 2019, 5:43 PM Stogner, Roy H <roy...@ic...>
wrote:

>
> On Tue, 1 Oct 2019, Renato Poli wrote:
>
> > I'm working with remeshing. Is there any helper class you I should look
> at?
> > I need to map a solution vector from a mesh to another, for exemple.
> > They share the same physical space, but the latter is much more refined
> in
> > a given area.
>
> See MeshFunction, combined with System::project_vector() or
> System::project_solution(), as in src/apps/projection.C
> ---
> Roy
>

Re: [Libmesh-users] Exchanging information between meshes

[Stogner, R. H <roy...@ic...>]

On Tue, 1 Oct 2019, Renato Poli wrote:

> I'm working with remeshing. Is there any helper class you I should look at?
> I need to map a solution vector from a mesh to another, for exemple.
> They share the same physical space, but the latter is much more refined in
> a given area.

See MeshFunction, combined with System::project_vector() or
System::project_solution(), as in src/apps/projection.C
---
Roy

Re: [Libmesh-users] Adapting systems_of_equations_ex7.C to heterogeneous Dirichlet BC

[Stogner, R. H <roy...@ic...>]

On Mon, 30 Sep 2019, Yuxiang Wang wrote:

> I am trying to adapt the systems_of_equations_ex7.C to heterogeneous
> Dirichlet BC, and I seem to hit a few challenges. Could you please help me
> out here?
>
> 1) My first attempt is simply adding the heterogeneous BC like a linear
> problem and used
> `dof_map.heterogenously_constrain_element_matrix_and_vector (Ke, Fe,
> dof_indices)` instead of the `dof_map.constrain_element_matrix_and_vector
> (Ke, Fe, dof_indices)` in the jacobian() function. Note that I did NOT
> have `system.rhs->add_vector (Fe, dof_indices)` in the jacobian() function.
> The solution seems to be slightly wrong - in the result file, some nodes
> did not have the right prescribed displacements, although their neighboring
> nodes seemed correct. Then I added the `system.rhs->add_vector (Fe,
> dof_indices)` in the jacobian() function then the results diverged.

Treating the heterogeneous BC like a linear problem is tricky if you
don't have control over what the solver's doing under the hood.

What we tend to do now in that case is a cheat: treat the constrained
DoFs as "mutable" in the C++ sense, do enforce_constraints_exactly()
before each residual or Jacobian evaluation, do homogeneous assembly
during the evaluations, then do enforce_constraints_exactly() one
final time after the solver returns.  See petsc_diff_solver.C for
example.

> 2) I noticed that the BCs have to be applied before the es.init().

> However, for NonlinearImplicitSystem we may want to divide the
> solution into multiple steps and we'd then need to update the
> heterogeneous BC per step.  Would this be possible?

If you have nonlinear (or time-varying) heterogeneous Dirichlet BCs,
System::reinit_constraints() can be used to recalculate the constraint
equations from them when you need.
---
Roy

[Libmesh-users] Exchanging information between meshes

[Renato P. <re...@gm...>]

Hi,

I'm working with remeshing. Is there any helper class you I should look at?
I need to map a solution vector from a mesh to another, for exemple.
They share the same physical space, but the latter is much more refined in
a given area.

Thanks.
Renato

[Libmesh-users] Adapting systems_of_equations_ex7.C to heterogeneous Dirichlet BC

[Yuxiang W. <yw...@vi...>]

Dear all,

I am trying to adapt the systems_of_equations_ex7.C to heterogeneous
Dirichlet BC, and I seem to hit a few challenges. Could you please help me
out here?

1) My first attempt is simply adding the heterogeneous BC like a linear
problem and used
`dof_map.heterogenously_constrain_element_matrix_and_vector (Ke, Fe,
dof_indices)` instead of the `dof_map.constrain_element_matrix_and_vector
(Ke, Fe, dof_indices)` in the jacobian() function. Note that I did NOT
have `system.rhs->add_vector (Fe, dof_indices)` in the jacobian() function.
The solution seems to be slightly wrong - in the result file, some nodes
did not have the right prescribed displacements, although their neighboring
nodes seemed correct. Then I added the `system.rhs->add_vector (Fe,
dof_indices)` in the jacobian() function then the results diverged.

2) I noticed that the BCs have to be applied before the es.init(). However,
for NonlinearImplicitSystem we may want to divide the solution into
multiple steps and we'd then need to update the heterogeneous BC per step.
Would this be possible?

Thank you!

Best,
Shawn

-- 
Yuxiang "Shawn" Wang, PhD
yw...@vi...
+1 (434) 284-0836

[Libmesh-users] evaluating solution field at arbitrary locations using subdivision

[Lee, J. H. <jae...@li...>]

Is there a way to evaluate solution field at arbitrary locations for subdivision FE? It seems like MeshFunction won’t work for subdivision FE because inverse_map is not implemented.
Also, it seems like the subdivision FE only supports the one-point Gaussian quadrature rule. Is this correct?
On a similar note, is there already a function in libMesh that will allow me to project solution field from one representation (subdivision) to another (lagrange)?
Thank you for your help!

Best,
Mike

Re: [Libmesh-users] Mesh refinement along a specific axis?

[Stogner, R. H <roy...@ic...>]

On Wed, 11 Sep 2019, Li Luo wrote:

> MeshFunction seems to manage the whole equation_systems.

It takes an EquationSystems reference, but you can query and/or
project only some of the solutions within.

> I wonder if there is any simpler operation that can transfer only
> the solutions between two different meshes? (with different
> partition, if possible). 

I'm afraid that's as simple as it gets; there's no anisotropic
uniformly_refine(), so in this case you'd have to actually generate a
new mesh.

> I wrote a routine that transfers the solutions between two nested
> meshes that one is uniformly refined from another (the refined mesh
> has the same partition as the coarse mesh). 

You know that you get this automatically in EquationSystems::reinit(),
right?
---
Roy

Re: [Libmesh-users] Mesh refinement along a specific axis?

[Li L. <li...@ka...>]

Thank you for your reply.

MeshFunction seems to manage the whole equation_systems. I wonder if there
is any simpler operation that can transfer only the solutions between two
different meshes? (with different partition, if possible).  And a cube
domain is enough for me.

I wrote a routine that transfers the solutions between two nested meshes
that one is uniformly refined from another (the refined mesh has the same
partition as the coarse mesh).
For some reason, I want the refined mesh to be denser than the coarse one
only along a particular direction, if libMesh provides another refining
operation similar to 'uniformly_refine()', then my code doesn't need to be
changed.

Best,
Li


On Wed, Sep 11, 2019 at 5:55 PM Stogner, Roy H <roy...@ic...>
wrote:

>
> On Wed, 11 Sep 2019, Li Luo wrote:
>
> > I used MeshTools::Generation::build_cube to generate a cartesian mesh.
> > I wonder is there any tool to refine the mesh along a specific axis
> > direction? For example, only refine uniformly along the z-axis.
>
> If you want to do it *uniformly*, on a cube, you can generate a new
> mesh with the increased resolution in that one direction, then use
> MeshFunction to project data from old mesh to new.
>
> If you want to do it adaptively, or on more general domains, you're
> out of luck.  This is something we've wanted to be able to do for years
> but we've got isotropic refinement assumptions too deeply held in too
> many places in the code.
> ---
> Roy
>


-- 

Postdoctoral Fellow
Extreme Computing Research Center
King Abdullah University of Science & Technology
https://sites.google.com/site/rolyliluo/

-- 

This message and its contents, including attachments are intended solely 
for the original recipient. If you are not the intended recipient or have 
received this message in error, please notify me immediately and delete 
this message from your computer system. Any unauthorized use or 
distribution is prohibited. Please consider the environment before printing 
this email.

Re: [Libmesh-users] Mesh refinement along a specific axis?

[Stogner, R. H <roy...@ic...>]

On Wed, 11 Sep 2019, Li Luo wrote:

> I used MeshTools::Generation::build_cube to generate a cartesian mesh.
> I wonder is there any tool to refine the mesh along a specific axis
> direction? For example, only refine uniformly along the z-axis.

If you want to do it *uniformly*, on a cube, you can generate a new
mesh with the increased resolution in that one direction, then use
MeshFunction to project data from old mesh to new.

If you want to do it adaptively, or on more general domains, you're
out of luck.  This is something we've wanted to be able to do for years
but we've got isotropic refinement assumptions too deeply held in too
many places in the code.
---
Roy

[Libmesh-users] Mesh refinement along a specific axis?

[Li L. <li...@ka...>]

Dear developer,

I used MeshTools::Generation::build_cube to generate a cartesian mesh.
I wonder is there any tool to refine the mesh along a specific axis
direction? For example, only refine uniformly along the z-axis.

Best,

Li Luo

-- 

This message and its contents, including attachments are intended solely 
for the original recipient. If you are not the intended recipient or have 
received this message in error, please notify me immediately and delete 
this message from your computer system. Any unauthorized use or 
distribution is prohibited. Please consider the environment before printing 
this email.

Re: [Libmesh-users] Mesh refinement and applying assemble function for weak form

[Stogner, R. H <roy...@ic...>]

On Tue, 10 Sep 2019, 서승진 wrote:

> Dear libmesh users,Hi. I currently use the mesh refinement class,
> and I have question about the mesh refinement.1) Before element is
> refined, it is active element,     but when the element is refined,
> it may not be active element, right?

Correct; once it's refined it becomes an "ancestor" element, with
children which are (until they're further refined or coarsened)
active.

> I currently use the assemble
> function, which approach to each element by using          " for
> (const auto & elem: mesh.active_local_element_ptr_range()) "  
> This "active_local_element_ptr_range()" may approach to the active
> element.         Is this loop can approach to the refined element
> (child element)?

After you've refined an element (and prepared the mesh for use again,
e.g. with EquationSystems::reinit()), that loop will now include the
child elements and will no longer include the original element.

> 2) Is it possible that specific elements which have refine
> flag is refined in several times?

If you manually flag elements, getting multiple levels of refinement
will require multiple "set flags, then reinit, then set more flags"
steps.  See our adaptivity or adjoints examples.
---
Roy

[Libmesh-users] Mesh refinement and applying assemble function for weak form

[<mis...@ka...>]

Dear libmesh users,Hi. I currently use the mesh refinement class, and I have question about the mesh refinement.1) Before element is refined, it is active element,     but when the element is refined, it may not be active element, right?    I currently use the assemble function, which approach to each element by using          " for (const auto & elem: mesh.active_local_element_ptr_range()) "   This "active_local_element_ptr_range()" may approach to the active element.         Is this loop can approach to the refined element (child element)?     Or, is it possible that the assemble function can set the weak form (matrix[Ke], vector[Fe]) to the refined elements?2) Is it possible that specific elements which have refine flag is refined in several times?         I am currently moving boundary node to simulate moving mesh,        and I should refine the boundary element for refinement of new generated area.    Maybe three times refinement of the boundary element can be sufficient for my case,     but when I conduct refinement in three times, there is no difference with one time refinement.    I can't find better option of mesh refinement class for this situation.Thanks for reading my problem!Best regard,Seungjin Seo

Re: [Libmesh-users] QTrap for TRI3SUBDIVISION

[Lee, J. H. <jae...@li...>]

I am having some trouble getting a stack trace, but I will continue to look more into what may the issue here. Another question — is there a function already available that support projecting solution from a system that uses subdivision onto a system that uses lagrange elements and vice versa?

On Aug 28, 2019, at 12:29 PM, John Peterson <jwp...@gm...<mailto:jwp...@gm...>> wrote:



On Wed, Aug 28, 2019 at 11:20 AM Lee, Jae Ho <jae...@li...<mailto:jae...@li...>> wrote:
John —

I just added TRI3SUBDIVISION to quadrature_trap_2D.C as a case, and I think it’s working now. I am still debugging things one step at a time, and one place that I got a complaint from is fe_lagrange.C — I may be confusing myself right now, but I think it’s complaining because TRI3SUBDIVISION is listed under FIRST order along with TRI3, but I think TRI3SUBDIVISION should be FOURTH order. I think I am confused because TRI3SUBDIVISION is listed in fe_lagrange.C.

Hmm... I guess somehow your code is calling lagrange_n_dofs_at_node() or lagrange_n_dofs(). That does seem wrong to me, unless somehow a LAGRANGE FE was created/reinit'd on a TRI3SUBDIVISION Elem. It might be helpful if you could give us a stack trace for the case where your code hits that line in fe_lagrange.C. Maybe we are building a Lagrange FE internally to the library for some reason...

Thanks for your report about the QTrap thing working, I will make a PR for that shortly.

--
John

Re: [Libmesh-users] QTrap for TRI3SUBDIVISION

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

On Wed, Aug 28, 2019 at 11:29 AM John Peterson <jwp...@gm...> wrote:

>
> Thanks for your report about the QTrap thing working, I will make a PR for
> that shortly.
>

https://github.com/libMesh/libmesh/pull/2230


-- 
John

Re: [Libmesh-users] QTrap for TRI3SUBDIVISION

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

On Wed, Aug 28, 2019 at 11:20 AM Lee, Jae Ho <jae...@li...> wrote:

> John —
>
> I just added TRI3SUBDIVISION to quadrature_trap_2D.C as a case, and I
> think it’s working now. I am still debugging things one step at a time, and
> one place that I got a complaint from is fe_lagrange.C — I may be confusing
> myself right now, but I think it’s complaining because TRI3SUBDIVISION is
> listed under FIRST order along with TRI3, but I think TRI3SUBDIVISION
> should be FOURTH order. I think I am confused because TRI3SUBDIVISION is
> listed in fe_lagrange.C.
>

Hmm... I guess somehow your code is calling lagrange_n_dofs_at_node()
or lagrange_n_dofs(). That does seem wrong to me, unless somehow a LAGRANGE
FE was created/reinit'd on a TRI3SUBDIVISION Elem. It might be helpful if
you could give us a stack trace for the case where your code hits that line
in fe_lagrange.C. Maybe we are building a Lagrange FE internally to the
library for some reason...

Thanks for your report about the QTrap thing working, I will make a PR for
that shortly.

-- 
John

Re: [Libmesh-users] QTrap for TRI3SUBDIVISION

[Lee, J. H. <jae...@li...>]

John —

I just added TRI3SUBDIVISION to quadrature_trap_2D.C as a case, and I think it’s working now. I am still debugging things one step at a time, and one place that I got a complaint from is fe_lagrange.C — I may be confusing myself right now, but I think it’s complaining because TRI3SUBDIVISION is listed under FIRST order along with TRI3, but I think TRI3SUBDIVISION should be FOURTH order. I think I am confused because TRI3SUBDIVISION is listed in fe_lagrange.C.

Thank you,
Mike

On Aug 28, 2019, at 12:05 PM, John Peterson <jwp...@gm...<mailto:jwp...@gm...>> wrote:

Hi Mike,

I think this (using QTrap with TRI3SUBDIVISION elements) should work just fine. Did you try it yet? Let us know and we can update the case statements for quadrature. There are many places like this in the library where TRI3 and TRI3SUBDIVISION elements can be used more or less interchangeably, so there's probably a bunch we haven't found yet.

--
John


On Wed, Aug 28, 2019 at 12:33 AM Lee, Jae Ho <jae...@li...<mailto:jae...@li...>> wrote:
Hello,

I am currently implementing something that makes use of TRI3SUBDIVISION elements, and one of the initial things that I want to do is to compute the nodal coordinates for each deformed element. I am trying to avoid the use of MeshFunction to do this. Of course since subdivision shape functions are not interpolatory, I can’t just use the solution of the coordinate system directly as the coordinates of the element vertices. So what I am trying to use QTrap that uses nodal quadrature points and evaluating the shape functions at the quadrature points for a given element to get something like

std::unique_ptr<libMesh::FEBase> fe(libMesh::FEBase::build(2, libMesh::SUBDIVISION));
std::unique_ptr<libMesh::QBase> qrule = libMesh::QBase::build(libMesh::QTRAP, 2, libMesh::FIRST);
fe->attach_quadrature_rule(qrule.get());
fe->reinit(elem);
for (unsigned int qp = 0; qp < qrule->n_points(); ++qp)
{
    for (unsigned int d = 0; d < 3; ++d)
    {
        for (unsigned int i = 0; i < phi.size(); ++i)
        {
            X_node[qp][d] += X_solution[i][d] * phi[i][qp];
        }
    }
}

X_solution is where the solution of the coordinate system has been stored.

Does this seem reasonable?
But quadrature_trap_2D.C complained because TRI3SUBDIVISION is not supported. Is it safe to hardcode TRI3SUBDIVISION as a case in the code? Or is there a different way to go about it (e.g. is there a different nodal quadrature rule that supports TRI3SUBDIVISION)?

Thank you,
Mike

Re: [Libmesh-users] QTrap for TRI3SUBDIVISION

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

Hi Mike,

I think this (using QTrap with TRI3SUBDIVISION elements) should work just
fine. Did you try it yet? Let us know and we can update the case statements
for quadrature. There are many places like this in the library where TRI3
and TRI3SUBDIVISION elements can be used more or less interchangeably, so
there's probably a bunch we haven't found yet.

--
John


On Wed, Aug 28, 2019 at 12:33 AM Lee, Jae Ho <jae...@li...> wrote:

> Hello,
>
> I am currently implementing something that makes use of TRI3SUBDIVISION
> elements, and one of the initial things that I want to do is to compute the
> nodal coordinates for each deformed element. I am trying to avoid the use
> of MeshFunction to do this. Of course since subdivision shape functions are
> not interpolatory, I can’t just use the solution of the coordinate system
> directly as the coordinates of the element vertices. So what I am trying to
> use QTrap that uses nodal quadrature points and evaluating the shape
> functions at the quadrature points for a given element to get something like
>
> std::unique_ptr<libMesh::FEBase> fe(libMesh::FEBase::build(2,
> libMesh::SUBDIVISION));
> std::unique_ptr<libMesh::QBase> qrule =
> libMesh::QBase::build(libMesh::QTRAP, 2, libMesh::FIRST);
> fe->attach_quadrature_rule(qrule.get());
> fe->reinit(elem);
> for (unsigned int qp = 0; qp < qrule->n_points(); ++qp)
> {
>     for (unsigned int d = 0; d < 3; ++d)
>     {
>         for (unsigned int i = 0; i < phi.size(); ++i)
>         {
>             X_node[qp][d] += X_solution[i][d] * phi[i][qp];
>         }
>     }
> }
>
> X_solution is where the solution of the coordinate system has been stored.
>
> Does this seem reasonable?
> But quadrature_trap_2D.C complained because TRI3SUBDIVISION is not
> supported. Is it safe to hardcode TRI3SUBDIVISION as a case in the code? Or
> is there a different way to go about it (e.g. is there a different nodal
> quadrature rule that supports TRI3SUBDIVISION)?
>
> Thank you,
> Mike
>
>

[Libmesh-users] QTrap for TRI3SUBDIVISION

[Lee, J. H. <jae...@li...>]

Hello,

I am currently implementing something that makes use of TRI3SUBDIVISION elements, and one of the initial things that I want to do is to compute the nodal coordinates for each deformed element. I am trying to avoid the use of MeshFunction to do this. Of course since subdivision shape functions are not interpolatory, I can’t just use the solution of the coordinate system directly as the coordinates of the element vertices. So what I am trying to use QTrap that uses nodal quadrature points and evaluating the shape functions at the quadrature points for a given element to get something like

std::unique_ptr<libMesh::FEBase> fe(libMesh::FEBase::build(2, libMesh::SUBDIVISION));
std::unique_ptr<libMesh::QBase> qrule = libMesh::QBase::build(libMesh::QTRAP, 2, libMesh::FIRST);
fe->attach_quadrature_rule(qrule.get());
fe->reinit(elem);
for (unsigned int qp = 0; qp < qrule->n_points(); ++qp)
{
    for (unsigned int d = 0; d < 3; ++d)
    {
        for (unsigned int i = 0; i < phi.size(); ++i)
        {
	    X_node[qp][d] += X_solution[i][d] * phi[i][qp];
        }
    }
}

X_solution is where the solution of the coordinate system has been stored.

Does this seem reasonable? 
But quadrature_trap_2D.C complained because TRI3SUBDIVISION is not supported. Is it safe to hardcode TRI3SUBDIVISION as a case in the code? Or is there a different way to go about it (e.g. is there a different nodal quadrature rule that supports TRI3SUBDIVISION)?

Thank you,
Mike

Re: [Libmesh-users] [RB EIM] Error estimator for EIM

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

Hello,

This email was empty, did you intend to send a question about RB EIM?

Best,
David

On Wed, Aug 21, 2019 at 5:32 AM <ss...@pu...> wrote:

>
> _______________________________________________
> Libmesh-users mailing list
> Lib...@li...
> https://lists.sourceforge.net/lists/listinfo/libmesh-users
>

[Libmesh-users] [RB EIM] Error estimator for EIM

[<ss...@pu...>]

Re: [Libmesh-users] Pure Co-located FVM-ALE Code: port to libMesh

[Stogner, R. H <roy...@ic...>]

On Thu, 15 Aug 2019, Andrew Parker wrote:

> I appreciate your help and comments, hope I didn’t cause any offense
> with the long list of questions.

No offense, but next time write to libmesh-users; you'll have more
eyes on them and the answers will be archived for future perusal.

I think John's pointer covered question 1 and to a lesser extent 4;
let me see if I can fill in the rest quickly.

>        1. I want to use the libMesh framework for the pure finite volume method, specifically, the cell-centered finite volume method with
>           all variables stored at cell centers, fluxes computed through faces, and face-connected nearest neighbours in the stencil in
>           serial and parallel (so level-1 one nearest face neighbours).  Do you perceive any problems in using libMesh for this?  Is there
>           an example to help me start?

> > I think you should take a close look at both the
> > "miscellaneous_ex5" and the new "vector_fe_ex5" examples, as those
> > both demonstrate using discontinuous Galerkin (DG) formulations,
> > which is going to be the best analog for implementing FVM in
> > libMesh. I think a lot of your questions will probably be answered
> > by looking through those examples, and further investigating some
> > of the APIs which they use. 

>        2. I mention ALE, because out-with the standard fixed grid cell-centered FVM method, mesh velocities stored at nodes will be
>           required.  Can I allocate, in parallel, and move in parallel, the mesh to allow for the ALE aspect of my work?  Again, any
>           examples?

We've had ALE user codes before, but there's no working support in the
library itself or examples; you'll have to handle mesh motion by hand.
We do have a couple sync functions if you just want to move local
nodes and let the library handle updating ghost node positions.

>        3. Starting with a mesh on disk I’d like to be able to do (for the above problem type), ghost cell creation, parallel partitioning,
>           graph reordering, mesh refinement and de-refinement, parallel comms, and parallel IO.  It appears libMesh has all of this, would
>           I be correct in assuming this?  From this point I’d want to define data on that parallel mesh, see next question.

Right.

>        4. I’d like to be able to allocate scalars, vectors, tensors in parallel on vertices, edges, faces, and cells (assuming 3D), or
>           vertices, edges and faces in 2D.  Is data allocation like that across the mesh possible? Do you have a simple example for say a
>           single cell in 2D/3D?  A particular question is with regards to what ghosting is done by libMesh for cells next to a processor
>           boundary (or face in 2D).  Is that “just” dealt with when requesting allocation of cell data?  For the nodes, faces, edges that
>           must be uniquely owned on a given processor, how does the allocation of data work on those entities when you want access to the
>           data allocated, on say, an edge on a processor that is not the owner?

We don't have independent allocation on vertices/edges/faces/cells;
each FE type defined which degrees of freedom go where.  If you're
trying to do something like DPG which defines some variables only on
interfaces between cells then you'll need to create a new FE type to
do so.

Our support for tensor-valued variables is only at the "bare shim"
level.  Unless you wanted to expand that, you'd need to store tensor
data as d^2 scalar variables (which is honestly also the best
supported way to store vector-valued variables, if your space is just
a tensor product of scalar-valued spaces and not an inherently
vector-valued FE type like Nedelec).

All degrees of freedom with support on local elements get ghosted; by
default so do DoFs supported on neighboring cells (which is useful for
e.g. gradient jump calculations) but you can narrow or widen that with
custom GhostingFunctor objects.

>        5. I’d like to be able to define boundary conditions on vertices, edges, and faces, and I’d like to be able to group cells.  I
>           would want to be able to do that from a mesh file (say MED/CGNS as I make extensive use of Salome) and have those groups
>           addressable (say as a collection of faces/cells) within the parallel mesh.  Is that possible?

>        6. Specifically, are subgroups of elements (nodes, edges, faces or cells) not just on the boundary but within the internals of the
>           mesh also a doable thing?

With cells, by subdomain id.  With faces, edges, and vertices, by
boundary id (which can be applied to "interior boundaries" despite the
oxymoron).

>        7. Are mesh metrics (cell centroid, cell volume, face centroid, face normal etc) computed by libMesh or me?  Including being
>           updated after refinement or mesh movement?

LibMesh gives methods for each of those.  They're computed on the fly,
not (except face normals from the FE object) cached within the library.

>        8. For some calculations I need to be able to refine only the edges (not the cells, in 2D), or the edges and then faces (but again
>           not the cells in 3D).  This produces hanging nodes and arbitrary polys, is this possible?

Not if you want us to handle mesh topology and constraint equations,
I'm afraid.  libMesh supports hanging nodes but only through isotropic
refinement, and doesn't support arbitrary polygons or polyhedra.

>        9. Does refinement do anything about the newly created data
>        (scalars/vectors) in the child elements (so refinement of an
>        edge with edge data, refinement of a face with face data,
>        refinement of a cell with cell data)? Same question on
>        de-refinement?

Child elements get data interpolated at nodes, then projected (holding
previously set data fixed) on edges, then faces, then interiors.

>        11. Outside of mesh velocities at nodes, I’d want to write
>        out cell-centered cell data, are those writers supported? 

Yes.

>        I would also want to write out face/edge-based solution data
>        on boundaries such as mass fluxes or shear stress, are those
>        writers supported?

To viz formats?  Only for simple interpretations, e.g. where you could
interpret boundary data as mid-edge+mid-face node values.

>        12. I make extensive use of CGNS, would this be a problem?

Not if you can write your own I/O classes.  ;-)

>        13. Is it possible to build a mesh from my own format?  So,
>        say I’ve read in the vertices, edges, faces, cells.  How
>        would I go about building a libMesh mesh, that is distributed
>        and has data allocated on the vertices, edges, faces, cells. 
>        This indeed may be a place to start as an example, in
>        particular being able to write back out that mesh in parallel
>        albeit for say two cells, in 2D, over two processors.  I’d
>        probably learn quite a bit from being able to do that.

Check out mesh_generation.C, or some of the small manually-constructed
meshes in our unit tests:

tests/base/overlapping_coupling_test.C
tests/mesh/boundary_info.C
tests/mesh/boundary_mesh.C
tests/mesh/boundary_points.C
tests/mesh/extra_integers.C
tests/mesh/mesh_function_dfem.C
tests/mesh/mixed_dim_mesh_test.C
tests/mesh/slit_mesh_test.C
tests/systems/equation_systems_test.C
tests/systems/systems_test.C
---
Roy