libmesh-users — General information for libmesh users

[Libmesh-users] Examples in topology optimization?

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

Dear all,

Sorry for the spam - I am curious in whether anyone knows about if there
are existing examples running topology optimization with libmesh. I googled
and found a few papers, but the source code seems to be unavailable.

Any hints would be appreciated!

Best,
Shawn

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

Re: [Libmesh-users] [RB] How to print output error bounds in the offline stage (RB training)

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

On Fri, Jan 18, 2019 at 1:37 AM <ss...@pu...> wrote:

> Hello, all.
>
>
>
> I have a question about printing output error bounds in the offline stage.
>
>
>
> We can see that maximum solution error bounds of each step in the RB
> training process, such as
>
> ==================================================================
>
> .
>
> ---- Basis dimension: 1 ----
>
> Performing RB solves on training set
>
> Maximum error bound is 1.07082
>
>
>
> Performing truth solve at parameter:
>
> .
>
> ==================================================================
>
>
>
> Besides, using the code "rb_eval.RB_output_error_bounds[0]" in the main
> program, I printed the output error bound after the RB training.
>
> However, I cannot find any way to print all output error bounds during an
> RB
> training.
>
> In other words, I want to see output error bound corresponding the basis
> dimension similar to the above code.
>
>
>
> Could you please tell me how to print them?
>
>
>
> I always thank you for your help.
>

I'm not really sure what you want to print... but the print statement you
referred to above is in rb_construction.C, where it does "libMesh::out <<
"Maximum error bound is " << training_greedy_error << std::endl <<
std::endl;". I guess you could do something like that in your code, or add
libMesh::out calls in your local libMesh in order to print more info?

Best,
David

[Libmesh-users] [RB] How to print output error bounds in the offline stage (RB training)

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

Hello, all.

 

I have a question about printing output error bounds in the offline stage.

 

We can see that maximum solution error bounds of each step in the RB
training process, such as

==================================================================

.

---- Basis dimension: 1 ----

Performing RB solves on training set

Maximum error bound is 1.07082

 

Performing truth solve at parameter:

.

==================================================================

 

Besides, using the code "rb_eval.RB_output_error_bounds[0]" in the main
program, I printed the output error bound after the RB training.

However, I cannot find any way to print all output error bounds during an RB
training.

In other words, I want to see output error bound corresponding the basis
dimension similar to the above code.

 

Could you please tell me how to print them?

 

I always thank you for your help.

 

Regards,

SKang

 

Re: [Libmesh-users] Question about the nodal solution value outputted via write_equation_systems() from L2_LAGRANGE elements

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

Hi David,

Ah got it - thanks a lot! I was attempting to
use write_discontinuous_equation_systems() and your help totally solved the
problem. Appreciate it!

Best,
Shawn

On Mon, Jan 14, 2019 at 6:46 PM David Knezevic <dav...@ak...>
wrote:

> On Mon, Jan 14, 2019 at 9:40 PM Yuxiang Wang <yw...@vi...> wrote:
>
>> Dear all,
>>
>> I am using the L2_LAGRANGE elements for stress recovery, because the
>> stress
>> are not continuous over the domain. I noticed that my
>> stress_system.solution->size() indicate that I have more than one value
>> for
>> each node and each stress component. This makes sense - I should have as
>> many values as the number of elements connecting to that node.
>>
>> However, when I used write_equation_system() functionality and visualize
>> in
>> ParaView, only one value is shown per node. I suspect that it is an
>> average
>> of all element DOFs on this node but am not sure - could you please help
>> share a hint?
>
>
> I'm not sure what paraview plots in this case, but you can plot a
> discontinuous solution by using write_discontinuous_exodusII, for example.
> If I recall correctly, that will create a new mesh with disconnected
> elements so that paraview will correctly plot it as discontinuous. The
> downside of this is the mesh no longer matches the mesh that you solved on
> (since the elements are disconnected) but that may be fine in your case.
>
> Best,
> David
>


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

Re: [Libmesh-users] Question about the nodal solution value outputted via write_equation_systems() from L2_LAGRANGE elements

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

On Mon, Jan 14, 2019 at 9:40 PM Yuxiang Wang <yw...@vi...> wrote:

> Dear all,
>
> I am using the L2_LAGRANGE elements for stress recovery, because the stress
> are not continuous over the domain. I noticed that my
> stress_system.solution->size() indicate that I have more than one value for
> each node and each stress component. This makes sense - I should have as
> many values as the number of elements connecting to that node.
>
> However, when I used write_equation_system() functionality and visualize in
> ParaView, only one value is shown per node. I suspect that it is an average
> of all element DOFs on this node but am not sure - could you please help
> share a hint?


I'm not sure what paraview plots in this case, but you can plot a
discontinuous solution by using write_discontinuous_exodusII, for example.
If I recall correctly, that will create a new mesh with disconnected
elements so that paraview will correctly plot it as discontinuous. The
downside of this is the mesh no longer matches the mesh that you solved on
(since the elements are disconnected) but that may be fine in your case.

Best,
David

[Libmesh-users] Question about the nodal solution value outputted via write_equation_systems() from L2_LAGRANGE elements

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

Dear all,

I am using the L2_LAGRANGE elements for stress recovery, because the stress
are not continuous over the domain. I noticed that my
stress_system.solution->size() indicate that I have more than one value for
each node and each stress component. This makes sense - I should have as
many values as the number of elements connecting to that node.

However, when I used write_equation_system() functionality and visualize in
ParaView, only one value is shown per node. I suspect that it is an average
of all element DOFs on this node but am not sure - could you please help
share a hint?

Best,
Shawn

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

Re: [Libmesh-users] Libmesh multi mesh for each rank / Libmesh moving boundary

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

Please use the libmesh-users list in the future; it's often responsive
even when an individual developer is busy.

On Tue, 8 Jan 2019, 서승진 wrote:

> Hi. I currently trying to solve the boundary deformation model in parallel.
> 
> 1) I have tried to write the code to solve the boundary deformation by using moving mesh,
> exactly saying, by moving the nodes in the boundary.
> 
> I used the MeshBase element iterator, and the boundary condition by using (elem->neighbor(side) == NULL) condition.
> Within the condition in the loop, I used
>   mesh.node(node_id) = Point (new_x, new_y, new_z);
> function, and it worked. However, the deformed boundary had spitz pattern. (Original shape of boundary was line)

I'm afraid I don't know what a "spitz pattern" is.  Could you upload a
picture somewhere?

> I want to make a model for deformed boundary without these spitz
> pattern, but it seems that I should link the moved boundary
> nodes each other by adding new boundary lines .

This seems very unlikely.  I might be misunderstanding the problem,
but generally with a deforming domain you want to change the geometry
but not the topology.

> 2) Libmesh can calculate the fem physics in one mesh in parallel, as
> like example introduction_ex4.C

> However, what I want to do is calculating two mesh with same fem
> physics in parallel.
> 
> By using Parallel::Communicator class with
> 
>   LibmeshInit init (argc,argv);
>   Parallel::Communicator & comm = init.comm();
> 
> I can assign the work for each processor along it's rank (ex, processor_id_type rank = comm.rank();)
> 
> Each processor will define mesh object, equation systems object, and system object, and then solve each system.

The key here is: how do you define the mesh object?  If you just use
comm above, then you will get a Mesh defined over all processors,
because init.comm() includes all processors libMesh sees.  You'll need
to do something like Communicator::split to get a new Communicator
defined over a subset of processors, then create a Mesh using just
that.

I don't think we have enough (or any?...) test coverage on that use
case, so please let us know ASAP (preferably via libmesh-devel or a
Github issue) if you run into anything which looks like a bug, or
(preferably via libmesh-users) if you need more assistance.
---
Roy

[Libmesh-users] New Update to XCode

[Caleb M P. <cal...@ut...>]

Hey all,

Came back from Christmas break and began running simulations and had an odd
"clang: error: linker command failed with exit code 1". Was sure that the
changes I made had no effect on this (simply changed a parameter value) and
then realized Xcode put out a new update 12/31. Last time I randomly had
this error, that was the issue. Am wondering if anyone else has these
issues every Xcode update for Mac?

Best,
Caleb

Re: [Libmesh-users] !libmesh_assert_pass' failed during an AMR operation

[Salazar De T. M. <sal...@ll...>]

On 1/4/19, 3:37 PM, "Stogner, Roy H" <roy...@ic...> wrote:

    
    On Thu, 3 Jan 2019, Salazar De Troya, Miguel via Libmesh-users wrote:
    
    >> Running an example with 10 processors (works in serial) fails with the error:
    
    >One of our examples?

    No. My own example so to speak. Sorry for the confusion.
    
    >> !libmesh_assert_pass' failed
    >>
    >> The program stops at libmesh_assert(test_unflagged(true)) in
    >> MeshRefinement::refine_and_coarsen_elements(). The remaining stack
    >> trace when TotalView stops is as follows (sorry for the rudimentary
    >> copy and paste, it’s all I can do now):
    >>
    >> MeshRefinement::test_unflagged (bool libmesh_dbg_var(libmesh_assert_pass))
    >> std::ostream & operator << (std::ostream & os, const Elem & e)
    >> Elem::print_info (std::ostream & os)
    >> std::string Elem::get_info ()
    >> bool DofObject::valid_id () <-- Stops here.
    
    >That's... baffling.  At that point _coarsen_elements should have
    >coarsened (and cleared flags of) anything flagged to coarsen, then
    >_refine_elements should have refined (and cleared flags of) anything
    >flagged to refine, then there should obviously be nothing left to
    >coarsen or to refine.
    
    >Hmm... the flag clearing in coarsening is actually done when the
    >parents get coarsen() called.  Is it possible that we're screwing up
    >make_coarsening_compatible somewhere and leaving a child marked
    >COARSEN without properly marking its parent COARSEN_INACTIVE?

   I noticed there are two test_unflagged() in MeshRefinement::refine_and_coarsen_elements(). It is the first one where the problem is (where the condition (coarsening_changed_mesh || refining_changed_mesh) is met).
    
    >> Please let me know if you need more details. I will keep investigating.
    
    >There's no output from that print_info?

   There is no more output. 

   I should apologize because the "stack trace" that I sent is not the one for the processor that failed (I'm a newbie at TotalView). This is the one for the processor that failed:

 libmesh_assert(!libmesh_assert_pass); in MeshRefinement::test_unflagged (bool libmesh_dbg_var(libmesh_assert_pass)) (mesh_refinement.C 467)
 libMesh::write_traceout() in MacroFunctions::report_error()
libMesh::print_trace(traceout);
gdb_worked = gdb_backtrace(out_stream);
int fd = mkstemp(temp_file);

and then I see two frames with what looks assembly language from the file libc.so.6

       
    >I don't know what I can do about this unless I can reproduce it.  I
    >could suggest some asserts that might catch the hypothesis above.

Could you please suggest those assertions? Sending the code could be a bit cumbersome.

Thanks
    ---
    Roy

Re: [Libmesh-users] !libmesh_assert_pass' failed during an AMR operation

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

On Thu, 3 Jan 2019, Salazar De Troya, Miguel via Libmesh-users wrote:

> Running an example with 10 processors (works in serial) fails with the error:

One of our examples?

> !libmesh_assert_pass' failed
>
> The program stops at libmesh_assert(test_unflagged(true)) in
> MeshRefinement::refine_and_coarsen_elements(). The remaining stack
> trace when TotalView stops is as follows (sorry for the rudimentary
> copy and paste, it’s all I can do now):
>
> MeshRefinement::test_unflagged (bool libmesh_dbg_var(libmesh_assert_pass))
> std::ostream & operator << (std::ostream & os, const Elem & e)
> Elem::print_info (std::ostream & os)
> std::string Elem::get_info ()
> bool DofObject::valid_id () <-- Stops here.

That's... baffling.  At that point _coarsen_elements should have
coarsened (and cleared flags of) anything flagged to coarsen, then
_refine_elements should have refined (and cleared flags of) anything
flagged to refine, then there should obviously be nothing left to
coarsen or to refine.

Hmm... the flag clearing in coarsening is actually done when the
parents get coarsen() called.  Is it possible that we're screwing up
make_coarsening_compatible somewhere and leaving a child marked
COARSEN without properly marking its parent COARSEN_INACTIVE?

> Please let me know if you need more details. I will keep investigating.

There's no output from that print_info?

I don't know what I can do about this unless I can reproduce it.  I
could suggest some asserts that might catch the hypothesis above.
---
Roy

[Libmesh-users] !libmesh_assert_pass' failed during an AMR operation

[Salazar De T. M. <sal...@ll...>]

Hello,

Running an example with 10 processors (works in serial) fails with the error:

!libmesh_assert_pass' failed

But I do not obtain any stack trace or more details. There are two files traceout_5_171629.txt and traceout_7_171631.txt which are empty. I run the code using TotalView and I get a bit more details. The program stops at libmesh_assert(test_unflagged(true)) in  MeshRefinement::refine_and_coarsen_elements(). The remaining stack trace when TotalView stops is as follows (sorry for the rudimentary copy and paste, it’s all I can do now):

MeshRefinement::test_unflagged (bool libmesh_dbg_var(libmesh_assert_pass))
std::ostream & operator << (std::ostream & os, const Elem & e)
Elem::print_info (std::ostream & os)
std::string Elem::get_info ()
bool DofObject::valid_id () <-- Stops here.

Please let me know if you need more details. I will keep investigating.

Thanks
Miguel

Miguel A. Salazar de Troya
Graduate Scholar, Lawrence Livermore National Laboratory
Ph.D. Candidate, University of Illinois at Urbana-Champaign
B141
Rm: 1085-5
Ph: 1(925) 422-6411

Re: [Libmesh-users] question about mesh related geometry problem

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

On Wed, 2 Jan 2019, Povolotskyi, Mykhailo wrote:

> Dear Libmesh developers and users,
> I need to solve the following problem.
> I have two subdomains in my mesh, e.g. 0 and 1.
> I have a set of quadrature points that belong to the subdomain 1.
> I need to find an approximate distance from each point to the subdomain 0.
> It will be enough if I can find an element in the subdomain 0 that is closest to a point in the subdomain 1.
> Do you have an idea how this can be done?

We don't have anything built into the library to help with that, I'm
afraid.  IIRC users in the past have used libMesh to construct distance
fields for use with a level-set method, though, and that's probably
what you want: construct the distance function (e.g. as a variable in
a separate ExplicitSystem on the same Mesh) with one fill sweep and
then reuse it over and over until/unless you modify the mesh.
---
Roy

[Libmesh-users] question about mesh related geometry problem

[Povolotskyi, M. <mpo...@pu...>]

Dear Libmesh developers and users,
I need to solve the following problem.
I have two subdomains in my mesh, e.g. 0 and 1.
I have a set of quadrature points that belong to the subdomain 1.
I need to find an approximate distance from each point to the subdomain 0.
It will be enough if I can find an element in the subdomain 0 that is closest to a point in the subdomain 1.
Do you have an idea how this can be done?
Thank you,
Michael.

Sent from Mail<https://go.microsoft.com/fwlink/?LinkId=550986> for Windows 10

Re: [Libmesh-users] output of a vector

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

On Thu, Dec 13, 2018 at 5:05 PM Michael Povolotskyi <mpo...@pu...>
wrote:

> Thank you!
>
> Do you mean that I have to store the system solution in some copy, then do
> what you suggested and restore the solution back?
>
If you need to keep the solution vector (e.g. for the next time step), then
yes, that's what I'd do. There may be another way though, maybe other can
weigh in on that.

David

Re: [Libmesh-users] output of a vector

[Michael P. <mpo...@pu...>]

Thank you!

Do you mean that I have to store the system solution in some copy, then 
do what you suggested and restore the solution back?


On 12/13/2018 04:54 PM, David Knezevic wrote:
> On Thu, Dec 13, 2018 at 4:35 PM Michael Povolotskyi 
> <mpo...@pu... <mailto:mpo...@pu...>> wrote:
>
>     Hello,
>
>     what is the recommended way to output a vector, previously added to a
>     system,  for visualization?
>
>     thank you!
>
>
> I normally do:
> *(system.solution) = system.get_vector("vec_name");
>
> and then write out the system as usual. This of course involves 
> copying the vector, I'm not sure if there's a way to avoid that copy 
> (usually the copy would be negligible compared to the rest of the 
> computations though...)
>
> David

Re: [Libmesh-users] output of a vector

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

On Thu, Dec 13, 2018 at 4:35 PM Michael Povolotskyi <mpo...@pu...>
wrote:

> Hello,
>
> what is the recommended way to output a vector, previously added to a
> system,  for visualization?
>
> thank you!


I normally do:
*(system.solution) = system.get_vector("vec_name");

and then write out the system as usual. This of course involves copying the
vector, I'm not sure if there's a way to avoid that copy (usually the copy
would be negligible compared to the rest of the computations though...)

David

[Libmesh-users] output of a vector

[Michael P. <mpo...@pu...>]

Hello,

what is the recommended way to output a vector, previously added to a 
system,  for visualization?

thank you!

Re: [Libmesh-users] Iterate all nodes in a certain BOUNDARY_ID

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

Hi Roy,

Awesome, thanks!

Best,
Shawn

On Wed, Dec 5, 2018 at 6:09 AM Stogner, Roy H <roy...@ic...>
wrote:

>
> On Tue, 4 Dec 2018, Yuxiang Wang wrote:
>
> > When I use `MeshTools::Generation::build_square`, I obtain a mesh with:
> >    enum BOUNDARY_ID {
> >      BOUNDARY_ID_MIN_Y = 0,
> >      BOUNDARY_ID_MAX_X = 1,
> >      BOUNDARY_ID_MAX_Y = 2,
> >      BOUNDARY_ID_MIN_X = 3,
> >    };
> >
> > However, I noticed that
> >
> >    MeshBase::const_node_iterator nodeit =
> > mesh.bid_nodes_begin(BOUNDARY_ID_MAX_X);
> >    const MeshBase::const_node_iterator node_end =
> > mesh.bid_nodes_end(BOUNDARY_ID_MAX_X);
> >
> > does not work for such IDs.
>
> Yes, that's correct.  build_square sets boundary ids on *sides*, but
> not on nodes.
>
> > Could you please help me in how to iterate through all the nodes?
>
> You could first use BoundaryInfo::build_node_list_from_side_list() to
> assign side boundary ids to the nodes on that side.  Corner nodes will
> be assigned 2 boundary ids.
> ---
> Roy
>


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

Re: [Libmesh-users] Elasticity example - single element result drastically different from ABAQUS

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

Ah OK - I found out that this has been a commonly asked question. In Abaqus
the C3D8 element "is base don Hughes' selective reduced integration (SRI)
method" and therefore not a really fully integrated element.

https://www.researchgate.net/post/Stiffness_matrix_of_abaqus_element_c3d8_and_by_my_8_noded_isoparametric_element_is_not_matching-can_anyone_help

Sorry again for the spam,
Shawn

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

On Tue, Dec 4, 2018, 22:59 Yuxiang Wang <yw...@vi... wrote:

> Sorry that I just realized that I wasn't clear in the boundary conditions:
>
> I fixed the MIN_X side of the model, and applied -Z traction on the
> surface of MAX_X. The resultant displacement on the MAX_X surface is about
> -3e-3 in ABAQUS, but only -2.2e-3 in the libmesh model. The other
> directions and the magnitude of the displacement vector is also drastically
> different.
>
> Shawn
>
> On Tue, Dec 4, 2018 at 10:53 PM Yuxiang Wang <yw...@vi...> wrote:
>
>> Dear all,
>>
>> Sorry to spam again.
>>
>> I was playing with ABAQUS & libmesh and noticed that the single element
>> result for simple linear elasticity model is different. I know that there
>> could be many factors that is leading to this (e.g. ABAQUS doing tricks
>> that I am not aware of), but I'd really appreciate it if anyone could help
>> provide some guesses on why this is the case, even for a simple
>> single-element model. Thanks!
>>
>> Both models are 2x2x2 cubes (from -1 to +1 so the coordinate matches the
>> bi-unit domain). I added a surface traction of -0.25 in the -Z direction,
>> and observe the resultant displacements in Z. The ABAQUS model is attached
>> (Job-1.inp), and the modified elasticity example is also attached. The
>> resultant displacement, especially at the Z direction, is totally different:
>>
>>
>> [image: image.png]
>>
>> Any tips would be much appreciated!
>>
>> Best,
>> Shawn
>> --
>> Yuxiang "Shawn" Wang, PhD
>> yw...@vi...
>> +1 (434) 284-0836
>>
>
>
> --
> Yuxiang "Shawn" Wang, PhD
> yw...@vi...
> +1 (434) 284-0836
>

Re: [Libmesh-users] Question about using NodeElems with GhostingFunctor

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

On Wed, Dec 5, 2018 at 8:51 AM Stogner, Roy H <roy...@ic...>
wrote:

>
> On Tue, 4 Dec 2018, David Knezevic wrote:
>
> > Yep, geometry is irrelevant for me. Also, I'm using ReplicatedMesh
> > currently, so I gather that geometric ghosting is not an issue
> > anyway.
>
> That's correct.
>
> > Hmm, ok... well this is interesting. I've used the approach I
> > described a lot on 3D meshes (e.g. for node-to-surface contact
> > couplings between surfaces of a 3D mesh) and I haven't run into any
> > problems with it so far. I do think I have tried cases where I have
> > NodeElems on domain boundaries where two processors' elements share
> > a Node, and I haven't had any issues with that.
>
> I guess it depends how you add the coupling terms.  If you're
> currently giving each NodeElem the same pid as its node, and if you're
> also exclusively *calculating* their interactions on that pid, then
> you should be fine.  I was thinking about the case where you specify
> interactions as integral terms along the element boundaries.
>

OK, that makes sense. I am exclusively calculating their interactions on
the same pid, so that explains why it works fine.



> > It would be interesting to reproduce the issue that you have in mind
> > in a test case to confirm that it really is a problem or not.
>
> It would.  4 quads in 2D should be sufficient.  This might be a
> good idea if you have time, even if it is currently working, just as
> defense against me accidentally breaking it in the future.  :-D
>

OK, I'll look into that at some stage, thanks for the suggestion.


> > The simplest change would be if I just coupled the 3D elements on
> > the surface directly rather than the nodes on the surface of those
> > elements.  That would be easy enough, but I figured that the
> > approach that I'm currently using is nicer since it doesn't
> > overallocated the sparsity pattern (e.g. there is no coupling needed
> > between the inner nodes on the elements on the surface).
>
> That's right.  If you were using integral terms (which is what you
> want for general contact problems, right, where the nodes might not
> line up perfectly?) then the boundary element approach should get the
> sparsity pattern right, though.
>

We currently do node-to-node (when the nodes line up perfectly) as well as
node-to-surface and surface-to-surface, which are two different approaches
that both handle the case where the nodes don't line up. Surface-to-surface
is the case you have in mind, and in that case I couple the 3D elements
directly (which is a bit of overallocation as discussed above, coupling
boundary elements as you described would be more precise). Node-to-surface
involves coupling a "slave node" with a "master side", and that's the case
where we use the NodeElem coupling, i.e. couple the slave NodeElem to the
master NodeElems.

Best,
David

Re: [Libmesh-users] Iterate all nodes in a certain BOUNDARY_ID

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

On Tue, 4 Dec 2018, Yuxiang Wang wrote:

> When I use `MeshTools::Generation::build_square`, I obtain a mesh with:
>    enum BOUNDARY_ID {
>      BOUNDARY_ID_MIN_Y = 0,
>      BOUNDARY_ID_MAX_X = 1,
>      BOUNDARY_ID_MAX_Y = 2,
>      BOUNDARY_ID_MIN_X = 3,
>    };
>
> However, I noticed that
>
>    MeshBase::const_node_iterator nodeit =
> mesh.bid_nodes_begin(BOUNDARY_ID_MAX_X);
>    const MeshBase::const_node_iterator node_end =
> mesh.bid_nodes_end(BOUNDARY_ID_MAX_X);
>
> does not work for such IDs.

Yes, that's correct.  build_square sets boundary ids on *sides*, but
not on nodes.

> Could you please help me in how to iterate through all the nodes?

You could first use BoundaryInfo::build_node_list_from_side_list() to
assign side boundary ids to the nodes on that side.  Corner nodes will
be assigned 2 boundary ids.
---
Roy

Re: [Libmesh-users] Question about using NodeElems with GhostingFunctor

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

On Tue, 4 Dec 2018, David Knezevic wrote:

> Yep, geometry is irrelevant for me. Also, I'm using ReplicatedMesh
> currently, so I gather that geometric ghosting is not an issue
> anyway.

That's correct.

> Hmm, ok... well this is interesting. I've used the approach I
> described a lot on 3D meshes (e.g. for node-to-surface contact
> couplings between surfaces of a 3D mesh) and I haven't run into any
> problems with it so far. I do think I have tried cases where I have
> NodeElems on domain boundaries where two processors' elements share
> a Node, and I haven't had any issues with that.

I guess it depends how you add the coupling terms.  If you're
currently giving each NodeElem the same pid as its node, and if you're
also exclusively *calculating* their interactions on that pid, then
you should be fine.  I was thinking about the case where you specify
interactions as integral terms along the element boundaries.

> It would be interesting to reproduce the issue that you have in mind
> in a test case to confirm that it really is a problem or not.

It would.  4 quads in 2D should be sufficient.  This might be a
good idea if you have time, even if it is currently working, just as
defense against me accidentally breaking it in the future.  :-D

> The simplest change would be if I just coupled the 3D elements on
> the surface directly rather than the nodes on the surface of those
> elements.  That would be easy enough, but I figured that the
> approach that I'm currently using is nicer since it doesn't
> overallocated the sparsity pattern (e.g. there is no coupling needed
> between the inner nodes on the elements on the surface).

That's right.  If you were using integral terms (which is what you
want for general contact problems, right, where the nodes might not
line up perfectly?) then the boundary element approach should get the
sparsity pattern right, though.
---
Roy

Re: [Libmesh-users] Elasticity example - single element result drastically different from ABAQUS

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

Sorry that I just realized that I wasn't clear in the boundary conditions:

I fixed the MIN_X side of the model, and applied -Z traction on the surface
of MAX_X. The resultant displacement on the MAX_X surface is about -3e-3 in
ABAQUS, but only -2.2e-3 in the libmesh model. The other directions and the
magnitude of the displacement vector is also drastically different.

Shawn

On Tue, Dec 4, 2018 at 10:53 PM Yuxiang Wang <yw...@vi...> wrote:

> Dear all,
>
> Sorry to spam again.
>
> I was playing with ABAQUS & libmesh and noticed that the single element
> result for simple linear elasticity model is different. I know that there
> could be many factors that is leading to this (e.g. ABAQUS doing tricks
> that I am not aware of), but I'd really appreciate it if anyone could help
> provide some guesses on why this is the case, even for a simple
> single-element model. Thanks!
>
> Both models are 2x2x2 cubes (from -1 to +1 so the coordinate matches the
> bi-unit domain). I added a surface traction of -0.25 in the -Z direction,
> and observe the resultant displacements in Z. The ABAQUS model is attached
> (Job-1.inp), and the modified elasticity example is also attached. The
> resultant displacement, especially at the Z direction, is totally different:
>
>
> [image: image.png]
>
> Any tips would be much appreciated!
>
> Best,
> Shawn
> --
> Yuxiang "Shawn" Wang, PhD
> yw...@vi...
> +1 (434) 284-0836
>


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

[Libmesh-users] Iterate all nodes in a certain BOUNDARY_ID

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

Dear all,

When I use `MeshTools::Generation::build_square`, I obtain a mesh with:
    enum BOUNDARY_ID {
      BOUNDARY_ID_MIN_Y = 0,
      BOUNDARY_ID_MAX_X = 1,
      BOUNDARY_ID_MAX_Y = 2,
      BOUNDARY_ID_MIN_X = 3,
    };

However, I noticed that

    MeshBase::const_node_iterator nodeit =
mesh.bid_nodes_begin(BOUNDARY_ID_MAX_X);
    const MeshBase::const_node_iterator node_end =
mesh.bid_nodes_end(BOUNDARY_ID_MAX_X);

does not work for such IDs. Could you please help me in how to iterate
through all the nodes?

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

Re: [Libmesh-users] Question about using NodeElems with GhostingFunctor

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

Hi Roy,

Thanks for your comments, I've replied below:

On Tue, Dec 4, 2018 at 5:01 PM Stogner, Roy H <roy...@ic...>
wrote:

>
> On Tue, 4 Dec 2018, David Knezevic wrote:
>
> > I use NodeElems with GhostingFunctor to couple unconnected nodes (i.e. I
> > attach NodeElems to the nodes that I want to couple), and this works
> well.
> > However, I realized that this approach requires me to set the
> processor_id
> > in the NodeElem to match the processor_id of the Node that it points to,
> as
> > per the code below:
> >
> >  for (const auto & elem : mesh.active_element_ptr_range())
> >  {
> >    if(elem->type() == NODEELEM)
> >    {
> >      elem->processor_id() = elem->node_ref(0).processor_id();
> >    }
> >  }
> >
> > This seems reasonable to me, but I just wanted to check if this is an
> > expected requirement?
>
> Huh.  I need to think about that.
>
> For a mesh to be consistent, each Node is required to have the same
> processor id as at least *one* of the elements which use it, but
> presumably you're tacking NodeElems onto an existing mesh which
> already satisfies that requirement, in which case the processor id of
> the NodeElems can be anything and still be valid.
>

Yes, that's right, I'm tacking NodeElems onto an existing mesh


> For geometric ghosting to work, the processor id of the NodeElem
> determines which processor can also see its ghosted elements.  But
> you're just using NodeElem as a hack to get just the right amount of
> coupling working, so geometry is irrelevant here I assume, in which
> case we still have no restrictions on NodeElem pid.
>

Yep, geometry is irrelevant for me. Also, I'm using ReplicatedMesh
currently, so I gather that geometric ghosting is not an issue anyway.



> For algebraic ghosting to work... we loop over active_local_elements
> when building the send_list... and that might be a problem for you,
> huh?!  It seems like not only would you need a NodeElem to have the
> right pid for that, but on domain boundaries where two processors'
> elements share a Node, you could need a NodeElem with *each* pid?!
>
> Well, either that or you'd need to use higher-dimensional boundary
> elements.  So then NodeElem would be the way to go for boundary
> coupling only on a 1-D mesh; you'd need Edge2/Edge3 elements on a 2D
> mesh, and quads or tris on a 3D mesh.  And each boundary element
> would need to match pid with its interior_parent() - it sounds like
> BoundaryInfo::add_elements() may be the safe way to work on these
> sorts of problems.
>

Hmm, ok... well this is interesting. I've used the approach I described a
lot on 3D meshes (e.g. for node-to-surface contact couplings between
surfaces of a 3D mesh) and I haven't run into any problems with it so far.
I do think I have tried cases where I have NodeElems on domain boundaries
where two processors'
elements share a Node, and I haven't had any issues with that. It would be
interesting to reproduce the issue that you have in mind in a test case to
confirm that it really is a problem or not.

The simplest change would be if I just coupled the 3D elements on the
surface directly rather than the nodes on the surface of those elements.
That would be easy enough, but I figured that the approach that I'm
currently using is nicer since it doesn't overallocated the sparsity
pattern (e.g. there is no coupling needed between the inner nodes on the
elements on the surface).

David