libmesh-users — General information for libmesh users

Re: [Libmesh-users] Integration points through the thickness for shell elements

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

Hi all,

To follow up my previous email, I think I can simplify my question as below:

It seems to me that if we have more than one integration point across the
shell thickness, the entire integration would need to be done by manually
constructing the Jacobian matrix in 3D to get the JxW values, without using
the 2D JxW entry in the current QUADSHELL4 element. This is because this
Jacobian in 2D does not take into account the nodal thicknesses and
directors, and there are no simple formula to use this information to
reconstruct the 3D JxW value. I tried to put these details in the attached
PDF in the previous email.

If you wouldn't mind, could you please share with me your thoughts on
whether my statement above is true? Or that I am missing something?

Again, I really appreciate your time to help me on this issue.

Best,
Shawn

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

On Fri, Nov 16, 2018, 22:56 Yuxiang Wang <yw...@vi... wrote:

> ( I am not sure whether the equations will display fine, so I also
> attached a PDF version of the text.)
>
>
>
> Hi all,
>
>
>
> Thank you all so much for the timely & quality responses! I really
> appreciate it. And my apologies for the delayed response.
>
>
>
> John – sorry that I did not explain myself clearly. I was using the
> notation where the mapping usually is from a geometry with thickness from
> physical domain  to the reference domain , although shape functions  and
> not a function of the third reference variable . As a result, I will
> still be integrating over all three axis.
>
>
>
> Ben & David - thank you so much again for the guidance. That said, I
> couldn’t figure out how to do the integration by adding a loop outside the
> existing quadrature point loop while utilizing existing JxW values (or more
> specifically, the determinant  part). I only know how to do that manually
> from scratch (i.e., calculate the Jacobian matrix, and manually calculate
> its determinant, times the quadrature weight). Would you kindly share
> whether you know a paper or textbook talking about this?
>
>
>
> My confusion, listed in detail, is below:
>
> For convenience I am using the notation in the Bathe textbook (which is
> consistent with his MITC4 original paper), where he used  to denote the
> shape functions of the node  (sometimes also referred to as  or  in other
> references), and  to denote the coordinates in the reference domain
> (sometimes also referred to as  in other references).
>
> Intuitively, I think I can do this:
>
> Where is the full 3D Jacobian determinant and the  is the surface
> Jacobian determinant that libmesh currently provides for QUADSHELL4, and  is
> the nodal thicknesses of node . However, I cannot mathematically prove
> this equation.
>
>
>
> For a general shell element, the coordinates of a point is given by
>
> The variable  denotes spatial dimensions of  and denotes the local node
> id. Then,  is the i-th component of the unit nodal director at node , and  is
> the nodal thicknesses of node  in the  direction. Note that the shape
> functions  only and is not a function of .
>
> Then, in general, the Jacobian matrix is ,
>
> And we can note the three columns as  and  (the curvilinear basis for the
> convected coordinate system). Also, .
>
> In the mean time, the Jacobian for the 2D element is
>
> And if we denote the two columns as  then if you don’t mind, I can write
> as .
>
>
>
> I went back to my intuition equation, the left hand side would be
>
> And the right hand side would be
>
> I couldn’t figure how whether in general . Intuitively, it seems to me
> that for variable thickness shells where , this equation does not hold.
>
>
>
> Thank you for your patience in reading through, and please let me know if
> my write-up is not clear enough. Any thoughts would be appreciated!
>
>
>
> Best,
>
> Shawn
>
>
>
> On Wed, Nov 14, 2018 at 8:52 AM David Knezevic <dav...@ak...>
> wrote:
>
>> On Wed, Nov 14, 2018 at 11:45 AM Benjamin W. Spencer via Libmesh-users <
>> lib...@li...> wrote:
>>
>>> It's pretty standard for shell elements to have multiple integration
>>> points through the thickness at every in-plane integration point.
>>> Integrating the response through the thickness allows you to represent the
>>> variation in the nonlinear constitutive response of the material through
>>> the cross-section, and come up with resultant quantities at the locations
>>> of the in-plane integration points, which are then integrated using
>>> standard procedures.
>>>
>>> I haven't really gotten too far into this yet, but I don't think
>>> accommodating those extra integration points would involve changing how the
>>> integration rules or data structures would work in libMesh. I think you
>>> would just evaluate vectors of properties at the standard integration
>>> points, with each entry in the vector representing a different point
>>> through the thickness. We are just getting started on the path of
>>> developing shell elements in  MOOSE, so our group will be looking into how
>>> to handle this.
>>>
>>
>> Yes, I agree with this description. We do this for modeling composites
>> for example, since they have different properties in each layer of the
>> composite which you can model via quadrature through the thickness.
>>
>> However, I think the libMesh example that is being discussed here is just
>> meant to be as simple as possible and hence it doesn't do this. Also, I
>> believe you can use analytical formulas for the integration through the
>> thickness in the case that the material is uniform, so I guess that is what
>> is done in the example.
>>
>> Best,
>> David
>>
>>
>>
>>> On 11/14/18, 8:19 AM, "John Peterson" <jwp...@gm...> wrote:
>>>
>>>     On Tue, Nov 13, 2018 at 11:22 PM Yuxiang Wang <yw...@vi...>
>>> wrote:
>>>
>>>     > Dear all,
>>>     >
>>>     > As one usually reads from literature (or commercial software
>>>     > documentation), usually, a shell element would need >= 2 Gaussian
>>>     > quadrature points through the thickness to capture its bending
>>> behavior.
>>>     > For example, in the LS-DYNA documentation
>>>     > <
>>> https://urldefense.proofpoint.com/v2/url?u=https-3A__www.dynasupport.com_tutorial_ls-2Ddyna-2Dusers-2Dguide_elements&d=DwICAg&c=54IZrppPQZKX9mLzcGdPfFD1hxrcB__aEkJFOKJFd00&r=hn5akMybrkn-1oiQB8nm_y7trT_BOQm9jBgbzQWwxXA&m=d16wjsgwuY6Xejdr47KKgE8srFi-kHjT92yv6KeNbt0&s=XxuqMBzy7V7pzqu5KNGEyqWuMwh-JQEAJerwVZsdQFU&e=>
>>> or
>>>     > mentioned in this paper
>>>     > <
>>>     >
>>> https://urldefense.proofpoint.com/v2/url?u=http-3A__web.mit.edu_kjb_www_Principal-5FPublications_Performance-5Fof-5Fthe-5FMITC3-2B-5Fand-5FMITC4-2B-5Fshell-5Felements-5Fin-5Fwidely-5Fused-5Fbenchmark-5Fproblems.pdf&d=DwICAg&c=54IZrppPQZKX9mLzcGdPfFD1hxrcB__aEkJFOKJFd00&r=hn5akMybrkn-1oiQB8nm_y7trT_BOQm9jBgbzQWwxXA&m=d16wjsgwuY6Xejdr47KKgE8srFi-kHjT92yv6KeNbt0&s=tSn1-9z_P8T0uME2jwoYDUO7AQEElPc8f3IjxytF-EA&e=
>>>     > >
>>>     > .
>>>     >
>>>
>>>     I  guess I'm confused about what you mean by "thickness". Our SHELL
>>>     elements are logically two-dimensional (have zero thickness) so IMO
>>> it
>>>     doesn't make sense ask about integration in the transverse
>>> direction...
>>>
>>>     --
>>>     John
>>>
>>>     _______________________________________________
>>>     Libmesh-users mailing list
>>>     Lib...@li...
>>>
>>> https://urldefense.proofpoint.com/v2/url?u=https-3A__lists.sourceforge.net_lists_listinfo_libmesh-2Dusers&d=DwICAg&c=54IZrppPQZKX9mLzcGdPfFD1hxrcB__aEkJFOKJFd00&r=hn5akMybrkn-1oiQB8nm_y7trT_BOQm9jBgbzQWwxXA&m=d16wjsgwuY6Xejdr47KKgE8srFi-kHjT92yv6KeNbt0&s=CrRFwn9nab2x8xz4ZBh1T-gH5-hRr1_hZ-mFsYImu-E&e=
>>>
>>>
>>>
>>> _______________________________________________
>>> Libmesh-users mailing list
>>> Lib...@li...
>>> https://lists.sourceforge.net/lists/listinfo/libmesh-users
>>>
>>
>
> --
> Yuxiang "Shawn" Wang, PhD
> yw...@vi...
> +1 (434) 284-0836
>

Re: [Libmesh-users] Integration points through the thickness for shell elements

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

( I am not sure whether the equations will display fine, so I also attached
a PDF version of the text.)



Hi all,



Thank you all so much for the timely & quality responses! I really
appreciate it. And my apologies for the delayed response.



John – sorry that I did not explain myself clearly. I was using the
notation where the mapping usually is from a geometry with thickness from
physical domain  to the reference domain , although shape functions  and
not a function of the third reference variable . As a result, I will still
be integrating over all three axis.



Ben & David - thank you so much again for the guidance. That said, I
couldn’t figure out how to do the integration by adding a loop outside the
existing quadrature point loop while utilizing existing JxW values (or more
specifically, the determinant  part). I only know how to do that manually
from scratch (i.e., calculate the Jacobian matrix, and manually calculate
its determinant, times the quadrature weight). Would you kindly share
whether you know a paper or textbook talking about this?



My confusion, listed in detail, is below:

For convenience I am using the notation in the Bathe textbook (which is
consistent with his MITC4 original paper), where he used  to denote the
shape functions of the node  (sometimes also referred to as  or  in other
references), and  to denote the coordinates in the reference domain
(sometimes also referred to as  in other references).

Intuitively, I think I can do this:

Where is the full 3D Jacobian determinant and the  is the surface Jacobian
determinant that libmesh currently provides for QUADSHELL4, and  is the
nodal thicknesses of node . However, I cannot mathematically prove this
equation.



For a general shell element, the coordinates of a point is given by

The variable  denotes spatial dimensions of  and denotes the local node id.
Then,  is the i-th component of the unit nodal director at node , and  is
the nodal thicknesses of node  in the  direction. Note that the shape
functions  only and is not a function of .

Then, in general, the Jacobian matrix is ,

And we can note the three columns as  and  (the curvilinear basis for the
convected coordinate system). Also, .

In the mean time, the Jacobian for the 2D element is

And if we denote the two columns as  then if you don’t mind, I can write as
.



I went back to my intuition equation, the left hand side would be

And the right hand side would be

I couldn’t figure how whether in general . Intuitively, it seems to me that
for variable thickness shells where , this equation does not hold.



Thank you for your patience in reading through, and please let me know if
my write-up is not clear enough. Any thoughts would be appreciated!



Best,

Shawn



On Wed, Nov 14, 2018 at 8:52 AM David Knezevic <dav...@ak...>
wrote:

> On Wed, Nov 14, 2018 at 11:45 AM Benjamin W. Spencer via Libmesh-users <
> lib...@li...> wrote:
>
>> It's pretty standard for shell elements to have multiple integration
>> points through the thickness at every in-plane integration point.
>> Integrating the response through the thickness allows you to represent the
>> variation in the nonlinear constitutive response of the material through
>> the cross-section, and come up with resultant quantities at the locations
>> of the in-plane integration points, which are then integrated using
>> standard procedures.
>>
>> I haven't really gotten too far into this yet, but I don't think
>> accommodating those extra integration points would involve changing how the
>> integration rules or data structures would work in libMesh. I think you
>> would just evaluate vectors of properties at the standard integration
>> points, with each entry in the vector representing a different point
>> through the thickness. We are just getting started on the path of
>> developing shell elements in  MOOSE, so our group will be looking into how
>> to handle this.
>>
>
> Yes, I agree with this description. We do this for modeling composites for
> example, since they have different properties in each layer of the
> composite which you can model via quadrature through the thickness.
>
> However, I think the libMesh example that is being discussed here is just
> meant to be as simple as possible and hence it doesn't do this. Also, I
> believe you can use analytical formulas for the integration through the
> thickness in the case that the material is uniform, so I guess that is what
> is done in the example.
>
> Best,
> David
>
>
>
>> On 11/14/18, 8:19 AM, "John Peterson" <jwp...@gm...> wrote:
>>
>>     On Tue, Nov 13, 2018 at 11:22 PM Yuxiang Wang <yw...@vi...>
>> wrote:
>>
>>     > Dear all,
>>     >
>>     > As one usually reads from literature (or commercial software
>>     > documentation), usually, a shell element would need >= 2 Gaussian
>>     > quadrature points through the thickness to capture its bending
>> behavior.
>>     > For example, in the LS-DYNA documentation
>>     > <
>> https://urldefense.proofpoint.com/v2/url?u=https-3A__www.dynasupport.com_tutorial_ls-2Ddyna-2Dusers-2Dguide_elements&d=DwICAg&c=54IZrppPQZKX9mLzcGdPfFD1hxrcB__aEkJFOKJFd00&r=hn5akMybrkn-1oiQB8nm_y7trT_BOQm9jBgbzQWwxXA&m=d16wjsgwuY6Xejdr47KKgE8srFi-kHjT92yv6KeNbt0&s=XxuqMBzy7V7pzqu5KNGEyqWuMwh-JQEAJerwVZsdQFU&e=>
>> or
>>     > mentioned in this paper
>>     > <
>>     >
>> https://urldefense.proofpoint.com/v2/url?u=http-3A__web.mit.edu_kjb_www_Principal-5FPublications_Performance-5Fof-5Fthe-5FMITC3-2B-5Fand-5FMITC4-2B-5Fshell-5Felements-5Fin-5Fwidely-5Fused-5Fbenchmark-5Fproblems.pdf&d=DwICAg&c=54IZrppPQZKX9mLzcGdPfFD1hxrcB__aEkJFOKJFd00&r=hn5akMybrkn-1oiQB8nm_y7trT_BOQm9jBgbzQWwxXA&m=d16wjsgwuY6Xejdr47KKgE8srFi-kHjT92yv6KeNbt0&s=tSn1-9z_P8T0uME2jwoYDUO7AQEElPc8f3IjxytF-EA&e=
>>     > >
>>     > .
>>     >
>>
>>     I  guess I'm confused about what you mean by "thickness". Our SHELL
>>     elements are logically two-dimensional (have zero thickness) so IMO it
>>     doesn't make sense ask about integration in the transverse
>> direction...
>>
>>     --
>>     John
>>
>>     _______________________________________________
>>     Libmesh-users mailing list
>>     Lib...@li...
>>
>> https://urldefense.proofpoint.com/v2/url?u=https-3A__lists.sourceforge.net_lists_listinfo_libmesh-2Dusers&d=DwICAg&c=54IZrppPQZKX9mLzcGdPfFD1hxrcB__aEkJFOKJFd00&r=hn5akMybrkn-1oiQB8nm_y7trT_BOQm9jBgbzQWwxXA&m=d16wjsgwuY6Xejdr47KKgE8srFi-kHjT92yv6KeNbt0&s=CrRFwn9nab2x8xz4ZBh1T-gH5-hRr1_hZ-mFsYImu-E&e=
>>
>>
>>
>> _______________________________________________
>> Libmesh-users mailing list
>> Lib...@li...
>> https://lists.sourceforge.net/lists/listinfo/libmesh-users
>>
>

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

[Libmesh-users] PetscDMWrapper functionality

[Miguel A. S. de T. <sal...@gm...>]

Hello,

How can I use this functionality? I tried simply passing the flag
--use_petsc_cm --node-major-dofs and the PETSc flags:

-ksp_type cg
-pc_type mg
-pc_mg_levels 3
-pc_mg_galerkin

but I am getting the error:

[0]PETSC ERROR: --------------------- Error Message
--------------------------------------------------------------
  1 [0]PETSC ERROR: *No support for this operation for this object type*
78  [0]PETSC ERROR: *This DM cannot coarsen*
  1 [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html
for trouble shooting.
  2 [0]PETSC ERROR: Petsc Release Version 3.8.3, Dec, 09, 2017
  3 [0]PETSC ERROR:
/g/g92/miguel/code/topsm/build-new/test/3D_linear_compliance/./3D_linear_compliance
on a petsc_fenics named quartz21 by miguel Fri Nov 16 12:16:02 2018
  4 [0]PETSC ERROR: Configure options --known-level1-dcache-size=32768
--known-level1-dcache-linesize=64 --known-level1-dcache-assoc=8
--known-sizeof-char=1 --known-sizeof-void-p=8 --known-sizeof-sh    ort=2
--known-sizeof-int=4 --known-sizeof-long=8 --known-sizeof-long-long=8
--known-sizeof-float=4 --known-sizeof-double=8 --known-sizeof-size_t=8
--known-bits-per-byte=8 --known-memcmp-ok=1 --k    nown-sizeof-MPI_Comm=4
--known-sizeof-MPI_Fint=4 --known-mpi-long-double=1 --known-mpi-int64_t=1
--known-mpi-c-double-complex=1 --known-sdot-returns-double=0
--known-snrm2-returns-double=0 --kno    wn-has-attribute-aligned=1
--download-hypre=yes --download-metis --download-ml=yes --download-mumps
--download-parmetis --download-ptscotch --download-scalapack
--download-suitesparse --download    -superilu
--known-mpi-shared-libraries=0 --with-batch --with-debugging=no
--with-language=cxx PETSC_ARCH=petsc_fenics --COPTFLAGS="-O3 -march=native
-mtune=native" --CXXOPTFLAGS="-O3 -march=nati    ve -mtune=native"
--FOPTFLAGS="-O3 -march=native -mtune=native"
  5 [0]PETSC ERROR: #1 DMCoarsen() line 2363 in
/usr/gapps/FEniCS/petsc-3.8.3/src/dm/interface/dm.c
  6 [0]PETSC ERROR: #2 PCSetUp_MG() line 618 in
/usr/gapps/FEniCS/petsc-3.8.3/src/ksp/pc/impls/mg/mg.c
  7 [0]PETSC ERROR: #3 PCSetUp() line 924 in
/usr/gapps/FEniCS/petsc-3.8.3/src/ksp/pc/interface/precon.c
  8 [0]PETSC ERROR: #4 KSPSetUp() line 381 in
/usr/gapps/FEniCS/petsc-3.8.3/src/ksp/ksp/interface/itfunc.c
  9 [0]PETSC ERROR: #5 KSPSolve() line 612 in
/usr/gapps/FEniCS/petsc-3.8.3/src/ksp/ksp/interface/itfunc.c
 10 [0]PETSC ERROR: #6 SNESSolve_NEWTONLS() line 224 in
/usr/gapps/FEniCS/petsc-3.8.3/src/snes/impls/ls/ls.c
 11 [0]PETSC ERROR: #7 SNESSolve() line 4109 in
/usr/gapps/FEniCS/petsc-3.8.3/src/snes/interface/snes.c

I do apply three levels of uniform refinement before solving my system.

Thanks
Miguel

Re: [Libmesh-users] Integration points through the thickness for shell elements

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

On Wed, Nov 14, 2018 at 11:45 AM Benjamin W. Spencer via Libmesh-users <
lib...@li...> wrote:

> It's pretty standard for shell elements to have multiple integration
> points through the thickness at every in-plane integration point.
> Integrating the response through the thickness allows you to represent the
> variation in the nonlinear constitutive response of the material through
> the cross-section, and come up with resultant quantities at the locations
> of the in-plane integration points, which are then integrated using
> standard procedures.
>
> I haven't really gotten too far into this yet, but I don't think
> accommodating those extra integration points would involve changing how the
> integration rules or data structures would work in libMesh. I think you
> would just evaluate vectors of properties at the standard integration
> points, with each entry in the vector representing a different point
> through the thickness. We are just getting started on the path of
> developing shell elements in  MOOSE, so our group will be looking into how
> to handle this.
>

Yes, I agree with this description. We do this for modeling composites for
example, since they have different properties in each layer of the
composite which you can model via quadrature through the thickness.

However, I think the libMesh example that is being discussed here is just
meant to be as simple as possible and hence it doesn't do this. Also, I
believe you can use analytical formulas for the integration through the
thickness in the case that the material is uniform, so I guess that is what
is done in the example.

Best,
David



> On 11/14/18, 8:19 AM, "John Peterson" <jwp...@gm...> wrote:
>
>     On Tue, Nov 13, 2018 at 11:22 PM Yuxiang Wang <yw...@vi...>
> wrote:
>
>     > Dear all,
>     >
>     > As one usually reads from literature (or commercial software
>     > documentation), usually, a shell element would need >= 2 Gaussian
>     > quadrature points through the thickness to capture its bending
> behavior.
>     > For example, in the LS-DYNA documentation
>     > <
> https://urldefense.proofpoint.com/v2/url?u=https-3A__www.dynasupport.com_tutorial_ls-2Ddyna-2Dusers-2Dguide_elements&d=DwICAg&c=54IZrppPQZKX9mLzcGdPfFD1hxrcB__aEkJFOKJFd00&r=hn5akMybrkn-1oiQB8nm_y7trT_BOQm9jBgbzQWwxXA&m=d16wjsgwuY6Xejdr47KKgE8srFi-kHjT92yv6KeNbt0&s=XxuqMBzy7V7pzqu5KNGEyqWuMwh-JQEAJerwVZsdQFU&e=>
> or
>     > mentioned in this paper
>     > <
>     >
> https://urldefense.proofpoint.com/v2/url?u=http-3A__web.mit.edu_kjb_www_Principal-5FPublications_Performance-5Fof-5Fthe-5FMITC3-2B-5Fand-5FMITC4-2B-5Fshell-5Felements-5Fin-5Fwidely-5Fused-5Fbenchmark-5Fproblems.pdf&d=DwICAg&c=54IZrppPQZKX9mLzcGdPfFD1hxrcB__aEkJFOKJFd00&r=hn5akMybrkn-1oiQB8nm_y7trT_BOQm9jBgbzQWwxXA&m=d16wjsgwuY6Xejdr47KKgE8srFi-kHjT92yv6KeNbt0&s=tSn1-9z_P8T0uME2jwoYDUO7AQEElPc8f3IjxytF-EA&e=
>     > >
>     > .
>     >
>
>     I  guess I'm confused about what you mean by "thickness". Our SHELL
>     elements are logically two-dimensional (have zero thickness) so IMO it
>     doesn't make sense ask about integration in the transverse direction...
>
>     --
>     John
>
>     _______________________________________________
>     Libmesh-users mailing list
>     Lib...@li...
>
> https://urldefense.proofpoint.com/v2/url?u=https-3A__lists.sourceforge.net_lists_listinfo_libmesh-2Dusers&d=DwICAg&c=54IZrppPQZKX9mLzcGdPfFD1hxrcB__aEkJFOKJFd00&r=hn5akMybrkn-1oiQB8nm_y7trT_BOQm9jBgbzQWwxXA&m=d16wjsgwuY6Xejdr47KKgE8srFi-kHjT92yv6KeNbt0&s=CrRFwn9nab2x8xz4ZBh1T-gH5-hRr1_hZ-mFsYImu-E&e=
>
>
>
> _______________________________________________
> Libmesh-users mailing list
> Lib...@li...
> https://lists.sourceforge.net/lists/listinfo/libmesh-users
>

Re: [Libmesh-users] Integration points through the thickness for shell elements

[Benjamin W. S. <ben...@in...>]

It's pretty standard for shell elements to have multiple integration points through the thickness at every in-plane integration point. Integrating the response through the thickness allows you to represent the variation in the nonlinear constitutive response of the material through the cross-section, and come up with resultant quantities at the locations of the in-plane integration points, which are then integrated using standard procedures.

I haven't really gotten too far into this yet, but I don't think accommodating those extra integration points would involve changing how the integration rules or data structures would work in libMesh. I think you would just evaluate vectors of properties at the standard integration points, with each entry in the vector representing a different point through the thickness. We are just getting started on the path of developing shell elements in  MOOSE, so our group will be looking into how to handle this.

-Ben

On 11/14/18, 8:19 AM, "John Peterson" <jwp...@gm...> wrote:

    On Tue, Nov 13, 2018 at 11:22 PM Yuxiang Wang <yw...@vi...> wrote:
    
    > Dear all,
    >
    > As one usually reads from literature (or commercial software
    > documentation), usually, a shell element would need >= 2 Gaussian
    > quadrature points through the thickness to capture its bending behavior.
    > For example, in the LS-DYNA documentation
    > <https://urldefense.proofpoint.com/v2/url?u=https-3A__www.dynasupport.com_tutorial_ls-2Ddyna-2Dusers-2Dguide_elements&d=DwICAg&c=54IZrppPQZKX9mLzcGdPfFD1hxrcB__aEkJFOKJFd00&r=hn5akMybrkn-1oiQB8nm_y7trT_BOQm9jBgbzQWwxXA&m=d16wjsgwuY6Xejdr47KKgE8srFi-kHjT92yv6KeNbt0&s=XxuqMBzy7V7pzqu5KNGEyqWuMwh-JQEAJerwVZsdQFU&e=> or
    > mentioned in this paper
    > <
    > https://urldefense.proofpoint.com/v2/url?u=http-3A__web.mit.edu_kjb_www_Principal-5FPublications_Performance-5Fof-5Fthe-5FMITC3-2B-5Fand-5FMITC4-2B-5Fshell-5Felements-5Fin-5Fwidely-5Fused-5Fbenchmark-5Fproblems.pdf&d=DwICAg&c=54IZrppPQZKX9mLzcGdPfFD1hxrcB__aEkJFOKJFd00&r=hn5akMybrkn-1oiQB8nm_y7trT_BOQm9jBgbzQWwxXA&m=d16wjsgwuY6Xejdr47KKgE8srFi-kHjT92yv6KeNbt0&s=tSn1-9z_P8T0uME2jwoYDUO7AQEElPc8f3IjxytF-EA&e=
    > >
    > .
    >
    
    I  guess I'm confused about what you mean by "thickness". Our SHELL
    elements are logically two-dimensional (have zero thickness) so IMO it
    doesn't make sense ask about integration in the transverse direction...
    
    -- 
    John
    
    _______________________________________________
    Libmesh-users mailing list
    Lib...@li...
    https://urldefense.proofpoint.com/v2/url?u=https-3A__lists.sourceforge.net_lists_listinfo_libmesh-2Dusers&d=DwICAg&c=54IZrppPQZKX9mLzcGdPfFD1hxrcB__aEkJFOKJFd00&r=hn5akMybrkn-1oiQB8nm_y7trT_BOQm9jBgbzQWwxXA&m=d16wjsgwuY6Xejdr47KKgE8srFi-kHjT92yv6KeNbt0&s=CrRFwn9nab2x8xz4ZBh1T-gH5-hRr1_hZ-mFsYImu-E&e=
    

Re: [Libmesh-users] Integration points through the thickness for shell elements

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

On Tue, Nov 13, 2018 at 11:22 PM Yuxiang Wang <yw...@vi...> wrote:

> Dear all,
>
> As one usually reads from literature (or commercial software
> documentation), usually, a shell element would need >= 2 Gaussian
> quadrature points through the thickness to capture its bending behavior.
> For example, in the LS-DYNA documentation
> <https://www.dynasupport.com/tutorial/ls-dyna-users-guide/elements> or
> mentioned in this paper
> <
> http://web.mit.edu/kjb/www/Principal_Publications/Performance_of_the_MITC3+_and_MITC4+_shell_elements_in_widely_used_benchmark_problems.pdf
> >
> .
>

I  guess I'm confused about what you mean by "thickness". Our SHELL
elements are logically two-dimensional (have zero thickness) so IMO it
doesn't make sense ask about integration in the transverse direction...

-- 
John

[Libmesh-users] Fwd: [Libmesh]: In parallel error :Segmentation Violation, probably memory access out of range

[Fatemeh C. S. <che...@gm...>]

Dear All,

I just created a public git repository for some section of the code which
causes the problem,

https://bitbucket.org/FatemehChe/reaction_diffusion/src/master/

Shortly speaking, this code is part of a PDE-Constrained optimization
problem, and this function is called at each iteration of an iterative
solver.
As I mentioned before it randomly gives us this error  "Segmentation
Violation, probably memory access out of range".
That's why I run it this function 10k times to show you the error, you can
see the error in the file called result/mono
Although, I am able to run it on my local machine (Mac/OS, CPU 2.5 GHz
Intel Core i7) with 4 cores, without any problem,
whenever I run the project on the cluster with more cores it gives me that
error after some iterations.

To run the project on the cluster,
- make
- sh benchmark

Thanks in advance.

Best regards,
Fatemeh

---------- Forwarded message ---------
From: John Peterson <jwp...@gm...>
Date: Tue, Nov 13, 2018 at 4:07 PM
Subject: Re: [Libmesh]: In parallel error :Segmentation Violation, probably
memory access out of range
To: Fatemeh Chegini Salzmann <che...@gm...>


Hi Fatemeh,

Two requests:

1.) please provide your entire application code if possible, we can't debug
code snippets
2.) please send help requests like this to lib...@li... or
open an issue on GitHub. That way you will be more likely to get help
people with more time and/or insight than me.

-- 
John


On Tue, Nov 13, 2018 at 3:38 AM Fatemeh Chegini Salzmann <
che...@gm...> wrote:

> Dear John,
>
> I have couple of transient PDEs  that need to be solved at each time
> iteration, where the solution of one PDE is used to solve the next PDE.
> As you can see shortly here:
> //
> -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
>     system_gateVariable.time  = 0;
>     system_diffusion.time        = 0;
>     system_reaction.time        = 0
>   unsigned int t_step;
>   for (t_step=1; t_step<(T); t_step++)
>   {
>
>     system_gateVariable.time += dt_T;
>     system_diffusion.time    += dt_T;
>     system_reaction.time     += dt_T;
>
>     *system_gateVariable.old_local_solution =
> *system_gateVariable.current_local_solution;
>     *system_reaction.old_local_solution =
> *system_reaction.current_local_solution;
>     *system_diffusion.old_local_solution =
> *system_diffusion.current_local_solution;
>
>     equation_systems.get_system("gateVariable").solve();
>     equation_systems.get_system("reaction").solve();
>
>   *system_diffusion.old_local_solution =
> *system_reaction.current_local_solution;
>    equation_systems.get_system("diffusion").solve();
> }
> //
> -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
>
>
>  and in the assembly function, I called the *old_solution* whenever i
> need to read the solution of each PDE. as follows:
> system.*old_solution*(dof_indices[l]);,
>
> I have this error when I run the project in parallel mode "*Caught signal
> number 11 SEGV: Segmentation Violation, probably memory access out of range*
> "
> I even tried to use *close()*  after solving each equation, e.g.*
> equation_systems.get_system("reaction").solution.close();* which didn't
> help.
> I don't know how to deal with this problem.
>
> I would appreciate any help/suggetion.
>
> Thank in advance,
> Fatemeh
>
> *FYI:*
>
> [1]PETSC ERROR: ------------------------------------------------
> ------------------------
>
> [1]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation,
> probably memory access out of range
>
> [1]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
>
> [1]PETSC ERROR: or see
> http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind
>
> [1]PETSC ERROR: or try http://valgrind.org on GNU/linux and Apple Mac OS
> X to find memory corruption errors
>
> [1]PETSC ERROR: configure using --with-debugging=yes, recompile, link, and
> run
>
> [1]PETSC ERROR: to get more information on the crash.
>
> [1]PETSC ERROR: --------------------- Error Message
> --------------------------------------------------------------
>
> [1]PETSC ERROR: Signal received
>
> [1]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html
> for trouble shooting.
>
> [1]PETSC ERROR: Petsc Release Version 3.8.3, Dec, 09, 2017
>
> [1]PETSC ERROR:
> /home/chegini/inverseProblem/LibmeshCode_InParallel_Cluster/./assemble on a
> arch-linux2-c-opt named icsnode17 by chegini Tue Nov 13 10:56:57 2018
>
> [1]PETSC ERROR: Configure options --prefix=/apps/petsc/3.8.3
> --download-hypre=1 --with-ssl=0 --with-debugging=no --with-pic=1
> --with-shared-libraries=1 --with-cc=mpicc --with-cxx=mpicxx
> --with-fc=mpif90 --download-fblaslapack=1 --download-metis=1
> --download-parmetis=1 --download-superlu_dist=1 --download-mumps=1
> --download-scalapack=1 --CC=mpicc --CXX=mpicxx --FC=mpif90 --F77=mpif77
> --F90=mpif90 --CFLAGS="-fPIC -fopenmp" --CXXFLAGS="-fPIC -fopenmp"
> --FFLAGS="-fPIC -fopenmp" --FCFLAGS="-fPIC -fopenmp" --F90FLAGS="-fPIC
> -fopenmp" --F77FLAGS="-fPIC -fopenmp"
> PETSC_DIR=/apps/petsc/3.8.3/src/petsc-3.8.3
>
> [1]PETSC ERROR: #1 User provided function() line 0 in  unknown file
>
> --------------------------------------------------------------------------
>
> MPI_ABORT was invoked on rank 1 in communicator MPI_COMM_WORLD
>
> with errorcode 59.
>
>
> NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
>
> You may or may not see output from other processes, depending on
>
> exactly when Open MPI kills them.
>
> --------------------------------------------------------------------------
>
> In: PMI_Abort(59, N/A)
>
> srun: Job step aborted: Waiting up to 32 seconds for job step to finish.
>
> slurmstepd: *** STEP 924595.0 ON icsnode17 CANCELLED AT
> 2018-11-13T11:00:55 ***
>
> srun: error: icsnode17: tasks 0,3,5-16,18-19: Killed
>
> srun: error: icsnode17: task 1: Exited with exit code 59
>
> srun: error: icsnode17: tasks 2,4,17: Killed
>

Re: [Libmesh-users] [Libmesh]: In parallel error :Segmentation Violation, probably memory access out of range

[Fatemeh C. S. <che...@gm...>]

Dear All,

I just created a public git repository for some section of the code which causes the problem,  

https://bitbucket.org/FatemehChe/reaction_diffusion/src/master/

Shortly speaking, this code is part of a PDE-Constrained optimization problem, and this function is called at each iteration of an iterative solver. 
As I mentioned before it randomly gives us this error  "Segmentation Violation, probably memory access out of range". 
That's why I run it this function 10k times to show you the error, you can see the error in the file called result/mono  
Although, I am able to run it on my local machine (Mac/OS, CPU 2.5 GHz Intel Core i7) with 4 cores, without any problem,  
whenever I run the project on the cluster with more cores it gives me that error after some iterations.

To run the project on the cluster,
- make
- sh benchmark

Thanks in advance.

Best regards,
Fatemeh



> 
> 
> ---------- Forwarded message ---------
> From: John Peterson <jwp...@gm... <mailto:jwp...@gm...>>
> Date: Tue, Nov 13, 2018 at 4:07 PM
> Subject: Re: [Libmesh]: In parallel error :Segmentation Violation, probably memory access out of range
> To: Fatemeh Chegini Salzmann <che...@gm... <mailto:che...@gm...>>
> 
> 
> Hi Fatemeh,
> 
> Two requests: 
> 
> 1.) please provide your entire application code if possible, we can't debug code snippets
> 2.) please send help requests like this to lib...@li... <mailto:lib...@li...> or open an issue on GitHub. That way you will be more likely to get help people with more time and/or insight than me.
> 
> -- 
> John
> 
> 
> On Tue, Nov 13, 2018 at 3:38 AM Fatemeh Chegini Salzmann <che...@gm... <mailto:che...@gm...>> wrote:
> Dear John,
> 
> I have couple of transient PDEs  that need to be solved at each time iteration, where the solution of one PDE is used to solve the next PDE.
> As you can see shortly here:
> // -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
>     system_gateVariable.time  = 0;
>     system_diffusion.time        = 0;  
>     system_reaction.time        = 0  
>   unsigned int t_step;
>   for (t_step=1; t_step<(T); t_step++)
>   {
> 
>     system_gateVariable.time += dt_T;
>     system_diffusion.time    += dt_T;  
>     system_reaction.time     += dt_T;  
>     
>     *system_gateVariable.old_local_solution = *system_gateVariable.current_local_solution; 
>     *system_reaction.old_local_solution = *system_reaction.current_local_solution;         
>     *system_diffusion.old_local_solution = *system_diffusion.current_local_solution;       
> 
>     equation_systems.get_system("gateVariable").solve(); 
>     equation_systems.get_system("reaction").solve(); 
>   
>   *system_diffusion.old_local_solution = *system_reaction.current_local_solution;       
>    equation_systems.get_system("diffusion").solve();
> }
> // -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------    
> 
>  and in the assembly function, I called the old_solution whenever i need to read the solution of each PDE. as follows:
> system.old_solution(dof_indices[l]);,   
> 
> I have this error when I run the project in parallel mode "Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range"
> I even tried to use close()  after solving each equation, e.g. equation_systems.get_system("reaction").solution.close(); which didn't help.
> I don't know how to deal with this problem.
> 
> I would appreciate any help/suggetion.
> 
> Thank in advance,
> Fatemeh
> 
> FYI:
> [1]PETSC ERROR: ------------------------------------------------------------------------
> [1]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range
> [1]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
> [1]PETSC ERROR: or see http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind <http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind>
> [1]PETSC ERROR: or try http://valgrind.org <http://valgrind.org/> on GNU/linux and Apple Mac OS X to find memory corruption errors
> [1]PETSC ERROR: configure using --with-debugging=yes, recompile, link, and run
> [1]PETSC ERROR: to get more information on the crash.
> [1]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
> [1]PETSC ERROR: Signal received
> [1]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html <http://www.mcs.anl.gov/petsc/documentation/faq.html> for trouble shooting.
> [1]PETSC ERROR: Petsc Release Version 3.8.3, Dec, 09, 2017
> [1]PETSC ERROR: /home/chegini/inverseProblem/LibmeshCode_InParallel_Cluster/./assemble on a arch-linux2-c-opt named icsnode17 by chegini Tue Nov 13 10:56:57 2018
> [1]PETSC ERROR: Configure options --prefix=/apps/petsc/3.8.3 --download-hypre=1 --with-ssl=0 --with-debugging=no --with-pic=1 --with-shared-libraries=1 --with-cc=mpicc --with-cxx=mpicxx --with-fc=mpif90 --download-fblaslapack=1 --download-metis=1 --download-parmetis=1 --download-superlu_dist=1 --download-mumps=1 --download-scalapack=1 --CC=mpicc --CXX=mpicxx --FC=mpif90 --F77=mpif77 --F90=mpif90 --CFLAGS="-fPIC -fopenmp" --CXXFLAGS="-fPIC -fopenmp" --FFLAGS="-fPIC -fopenmp" --FCFLAGS="-fPIC -fopenmp" --F90FLAGS="-fPIC -fopenmp" --F77FLAGS="-fPIC -fopenmp" PETSC_DIR=/apps/petsc/3.8.3/src/petsc-3.8.3
> [1]PETSC ERROR: #1 User provided function() line 0 in  unknown file
> --------------------------------------------------------------------------
> MPI_ABORT was invoked on rank 1 in communicator MPI_COMM_WORLD
> with errorcode 59.
> 
> NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
> You may or may not see output from other processes, depending on
> exactly when Open MPI kills them.
> --------------------------------------------------------------------------
> In: PMI_Abort(59, N/A)
> srun: Job step aborted: Waiting up to 32 seconds for job step to finish.
> slurmstepd: *** STEP 924595.0 ON icsnode17 CANCELLED AT 2018-11-13T11:00:55 ***
> srun: error: icsnode17: tasks 0,3,5-16,18-19: Killed
> srun: error: icsnode17: task 1: Exited with exit code 59
> srun: error: icsnode17: tasks 2,4,17: Killed
> 
> 
> 
> 

Re: [Libmesh-users] Integration points through the thickness for shell elements

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

Correction: "small-strain assumptions " should be "thin-shell assumptions".
Sorry.



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

On Tue, Nov 13, 2018, 22:21 Yuxiang Wang <yw...@vi... wrote:

> Dear all,
>
> As one usually reads from literature (or commercial software
> documentation), usually, a shell element would need >= 2 Gaussian
> quadrature points through the thickness to capture its bending behavior.
> For example, in the LS-DYNA documentation
> <https://www.dynasupport.com/tutorial/ls-dyna-users-guide/elements> or
> mentioned in this paper
> <http://web.mit.edu/kjb/www/Principal_Publications/Performance_of_the_MITC3+_and_MITC4+_shell_elements_in_widely_used_benchmark_problems.pdf>
> .
>
> I noted that the QUADSHELL4 and QUADSHELL8 elements have only the 4
> quadrature point at the zeta=0 plane. I was thinking about manually adding
> one more layer for numerical integration, but just wonder that - would it
> make sense to build that in QUADSHELL4 or QUADSHELL8 default qrule? Or
> would you think it would be cleaner to manually add one more loop?
>
> PS: The existing MITC4 shell example (actually Q4Gamma24 element) seems to
> be using only one layer of quadrature points across the thickness, and the
> results does not match ABAQUS well with thicker shells. I am still trying
> to understand whether it is because of some small-strain assumptions during
> implementation of the Q4Gamma24 theory, or whether it is because of the
> 1-layer integration points through the thickness. I'll report back when I
> have an idea!
>
> PPS: Reading a French textbook with Google Translate is really
> challenging... :)
>
> Best,
> Shawn
>
> --
> Yuxiang "Shawn" Wang, PhD
> yw...@vi...
> +1 (434) 284-0836
>

[Libmesh-users] Integration points through the thickness for shell elements

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

Dear all,

As one usually reads from literature (or commercial software
documentation), usually, a shell element would need >= 2 Gaussian
quadrature points through the thickness to capture its bending behavior.
For example, in the LS-DYNA documentation
<https://www.dynasupport.com/tutorial/ls-dyna-users-guide/elements> or
mentioned in this paper
<http://web.mit.edu/kjb/www/Principal_Publications/Performance_of_the_MITC3+_and_MITC4+_shell_elements_in_widely_used_benchmark_problems.pdf>
.

I noted that the QUADSHELL4 and QUADSHELL8 elements have only the 4
quadrature point at the zeta=0 plane. I was thinking about manually adding
one more layer for numerical integration, but just wonder that - would it
make sense to build that in QUADSHELL4 or QUADSHELL8 default qrule? Or
would you think it would be cleaner to manually add one more loop?

PS: The existing MITC4 shell example (actually Q4Gamma24 element) seems to
be using only one layer of quadrature points across the thickness, and the
results does not match ABAQUS well with thicker shells. I am still trying
to understand whether it is because of some small-strain assumptions during
implementation of the Q4Gamma24 theory, or whether it is because of the
1-layer integration points through the thickness. I'll report back when I
have an idea!

PPS: Reading a French textbook with Google Translate is really
challenging... :)

Best,
Shawn

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

Re: [Libmesh-users] Get the n_nodes for a certain bid?

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

Thank you Roy & John! That's great help. Appreciate it :)

Best,
Shawn

On Mon, Nov 12, 2018 at 4:21 PM John Peterson <jwp...@gm...> wrote:

>
>
> On Mon, Nov 12, 2018 at 4:31 PM Stogner, Roy H <roy...@ic...>
> wrote:
>
>>
>> On Mon, 12 Nov 2018, Yuxiang Wang wrote:
>>
>> > Quick question: how do I get the n_nodes for a certain boundary set
>> (given
>> > bid)? I know I can do:
>> >
>> > MeshBase::const_node_iterator nodeit =
>> > mesh.bid_nodes_begin(NODE_BOUNDARY_ID);
>> > const MeshBase::const_node_iterator node_end =
>> > mesh.bid_nodes_end(NODE_BOUNDARY_ID);
>> >
>> > to loop through all the nodes, but would there be a way to find out how
>> > many nodes are there in the same set?
>>
>> Just std::distance, I believe.
>>
>
> This might be a bit slow, depending on what else you are doing. Another
> possibility would be to call:
>
> std::vector<std::tuple<dof_id_type, boundary_id_type>> vec
>   = mesh.get_boundary_info().build_node_list();
>
> which returns a vector of _all_ the (node_id, boundary_id) pairs in the
> mesh (they will be sorted by node_id). You could then count up the number
> of entries with a specific boundary_id as needed...
>
> --
> John
>


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

Re: [Libmesh-users] Get the n_nodes for a certain bid?

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

On Mon, Nov 12, 2018 at 4:31 PM Stogner, Roy H <roy...@ic...>
wrote:

>
> On Mon, 12 Nov 2018, Yuxiang Wang wrote:
>
> > Quick question: how do I get the n_nodes for a certain boundary set
> (given
> > bid)? I know I can do:
> >
> > MeshBase::const_node_iterator nodeit =
> > mesh.bid_nodes_begin(NODE_BOUNDARY_ID);
> > const MeshBase::const_node_iterator node_end =
> > mesh.bid_nodes_end(NODE_BOUNDARY_ID);
> >
> > to loop through all the nodes, but would there be a way to find out how
> > many nodes are there in the same set?
>
> Just std::distance, I believe.
>

This might be a bit slow, depending on what else you are doing. Another
possibility would be to call:

std::vector<std::tuple<dof_id_type, boundary_id_type>> vec
  = mesh.get_boundary_info().build_node_list();

which returns a vector of _all_ the (node_id, boundary_id) pairs in the
mesh (they will be sorted by node_id). You could then count up the number
of entries with a specific boundary_id as needed...

-- 
John

Re: [Libmesh-users] Get the n_nodes for a certain bid?

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

On Mon, 12 Nov 2018, Yuxiang Wang wrote:

> Quick question: how do I get the n_nodes for a certain boundary set (given
> bid)? I know I can do:
>
> MeshBase::const_node_iterator nodeit =
> mesh.bid_nodes_begin(NODE_BOUNDARY_ID);
> const MeshBase::const_node_iterator node_end =
> mesh.bid_nodes_end(NODE_BOUNDARY_ID);
>
> to loop through all the nodes, but would there be a way to find out how
> many nodes are there in the same set?

Just std::distance, I believe.
---
Roy

[Libmesh-users] Get the n_nodes for a certain bid?

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

Dear all,

Sorry for the spam again.

Quick question: how do I get the n_nodes for a certain boundary set (given
bid)? I know I can do:

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

to loop through all the nodes, but would there be a way to find out how
many nodes are there in the same set?

Best,
Shawn

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

Re: [Libmesh-users] A node does not have dof attached to it?

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

Ah - thank you for the hint!

Totally my bad. I accidentally put L2_LAGRANGE while the element should
really be LAGRANGE. That makes perfect sense now.

I have attached my main function just in case someone else in the future
may find this reference.

Thank you so much Roy!

Best,
Shawn

// Begin the main program.
int main (int argc, char ** argv) {
// Initialize libMesh.
LibMeshInit init (argc, argv);

Mesh mesh(init.comm(), 3);

MeshTools::Generation::build_square (mesh,
2, 2,
-1., 1.,
-1., 1.,
QUAD4);

mesh.set_spatial_dimension (3);

EquationSystems equation_systems (mesh);

ExplicitSystem & director_system = equation_systems.add_system<ExplicitSystem>
("DirectorSystem");
director_system.add_variable("Vnx", FIRST, LAGRANGE);
director_system.add_variable("Vny", FIRST, L2_LAGRANGE);

equation_systems.init ();


// Use an explicit system to store nodal directors
const unsigned int director_system_id = director_system.number();


for (auto node : mesh.node_ptr_range()) {
std::cout << node->n_comp(director_system_id, 0) << std::endl;
std::cout << node->n_comp(director_system_id, 1) << std::endl;
}

return 0;
}

And the output:

./example-dbg
1
0
1
0
1
0
1
0
1
0
1
0
1
0
1
0
1
0


On Thu, Nov 8, 2018 at 6:13 AM Stogner, Roy H <roy...@ic...>
wrote:

>
> On Thu, 8 Nov 2018, Yuxiang Wang wrote:
>
> > Sorry for the naive question. In short, I have a system with three
> > variables (first order Lagrangian elements). I thought that for each of
> the
> > node in the mesh, it would have a DOF attached to it with n_comp==1.
> > However, when I used the following code snippet to check, the n_comp is
> > actually 0 (not 1).
> >
> > for (auto node : mesh.node_ptr_range())
> >  std::cout << node->n_comp(0, 0) << std::endl;
> >
> > What I was trying to achieve is to iterate through all the nodes, and
> > assign solution values to each variable DOF attached to the node from a
> > Eigen::Vector3d. Then I got the error message about:
> >
> > Assertion `comp < this->n_comp(s,var)' failed.
> > comp = 0
> > this->n_comp(s,var) = 0
> >
> > Any help would be appreciated. Thank you!
>
> Are you on a mesh with only first-order geometric elements?  If not
> then the higher-order nodes (e.g. the mid-edge nodes on a TRI6) will
> have n_comp == 0 for first-order Lagrange variables.
> ---
> Roy
>


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

Re: [Libmesh-users] A node does not have dof attached to it?

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

On Thu, 8 Nov 2018, Yuxiang Wang wrote:

> Sorry for the naive question. In short, I have a system with three
> variables (first order Lagrangian elements). I thought that for each of the
> node in the mesh, it would have a DOF attached to it with n_comp==1.
> However, when I used the following code snippet to check, the n_comp is
> actually 0 (not 1).
>
> for (auto node : mesh.node_ptr_range())
>  std::cout << node->n_comp(0, 0) << std::endl;
>
> What I was trying to achieve is to iterate through all the nodes, and
> assign solution values to each variable DOF attached to the node from a
> Eigen::Vector3d. Then I got the error message about:
>
> Assertion `comp < this->n_comp(s,var)' failed.
> comp = 0
> this->n_comp(s,var) = 0
>
> Any help would be appreciated. Thank you!

Are you on a mesh with only first-order geometric elements?  If not
then the higher-order nodes (e.g. the mid-edge nodes on a TRI6) will
have n_comp == 0 for first-order Lagrange variables.
---
Roy

[Libmesh-users] A node does not have dof attached to it?

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

Dear all,

Sorry for the naive question. In short, I have a system with three
variables (first order Lagrangian elements). I thought that for each of the
node in the mesh, it would have a DOF attached to it with n_comp==1.
However, when I used the following code snippet to check, the n_comp is
actually 0 (not 1).

for (auto node : mesh.node_ptr_range())
  std::cout << node->n_comp(0, 0) << std::endl;

What I was trying to achieve is to iterate through all the nodes, and
assign solution values to each variable DOF attached to the node from a
Eigen::Vector3d. Then I got the error message about:

Assertion `comp < this->n_comp(s,var)' failed.
comp = 0
this->n_comp(s,var) = 0

Any help would be appreciated. Thank you!

Best,
Shawn

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

Re: [Libmesh-users] question on Trapezoidal integration rule

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

On 11/07/2018 04:02 PM, Stogner, Roy H wrote:
> On Wed, 7 Nov 2018, Michael Povolotskyi wrote:
>
>> On 11/07/2018 02:59 PM, Stogner, Roy H wrote:
>>>   On Wed, 7 Nov 2018, Michael Povolotskyi wrote:
>>>
>>>>   if I use the trapezoidal rule:
>>>>
>>>>   QTrap qrule(dim)
>>>>
>>>>   Can I assume the quadrature points are ordered in the same way at the
>>>>   element
>>>>   nodes?
>>>   Doesn't look that way, in general, at first glance.  We build the quad
>>>   and hex QTrap rules via tensor products, but that's not how we order
>>>   nodes.
>>>   ---
>>>   Roy
>> Thank you,
>> then how to solve the following problem.
>>
>> I have an element that has n_nodes;
>> It has the same number of quadrature points that are at the nodal positions
>> (trapezoidal rule is assumed).
>> How to find the node number that corresponds to a quadrature point?
>> Should I compare coordinates?
> If you do that at runtime then it will throw a ton of unnecessary
> calculation into an inner loop.  You could just hardcode the
> permutation vectors; looking over quadrature_trap*.C it looks like
> (mapping from quadrature point numbers to node numbers, and don't
> trust me on this without verifying it):
>
> QUAD4: {0, 1, 3, 2}
> HEX8:  {0, 1, 3, 2, 4, 5, 7, 6}
>
> And then edge/triangle/prism elements are all the identity
> permutation... and we don't support QTrap on pyramids apparently?
> ---
> Roy
Thank you,
I'll code a permutation vector inside my application.
Michael.

Re: [Libmesh-users] question on Trapezoidal integration rule

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

On Wed, 7 Nov 2018, Michael Povolotskyi wrote:

> On 11/07/2018 02:59 PM, Stogner, Roy H wrote:
>>  On Wed, 7 Nov 2018, Michael Povolotskyi wrote:
>>
>>>  if I use the trapezoidal rule:
>>>
>>>  QTrap qrule(dim)
>>>
>>>  Can I assume the quadrature points are ordered in the same way at the
>>>  element
>>>  nodes?
>>  Doesn't look that way, in general, at first glance.  We build the quad
>>  and hex QTrap rules via tensor products, but that's not how we order
>>  nodes.
>>  ---
>>  Roy
> Thank you,
> then how to solve the following problem.
>
> I have an element that has n_nodes;
> It has the same number of quadrature points that are at the nodal positions 
> (trapezoidal rule is assumed).
> How to find the node number that corresponds to a quadrature point?
> Should I compare coordinates?

If you do that at runtime then it will throw a ton of unnecessary
calculation into an inner loop.  You could just hardcode the
permutation vectors; looking over quadrature_trap*.C it looks like
(mapping from quadrature point numbers to node numbers, and don't
trust me on this without verifying it):

QUAD4: {0, 1, 3, 2}
HEX8:  {0, 1, 3, 2, 4, 5, 7, 6}

And then edge/triangle/prism elements are all the identity
permutation... and we don't support QTrap on pyramids apparently?
---
Roy

Re: [Libmesh-users] question on Trapezoidal integration rule

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

On 11/07/2018 02:59 PM, Stogner, Roy H wrote:
> On Wed, 7 Nov 2018, Michael Povolotskyi wrote:
>
>> if I use the trapezoidal rule:
>>
>> QTrap qrule(dim)
>>
>> Can I assume the quadrature points are ordered in the same way at the element
>> nodes?
> Doesn't look that way, in general, at first glance.  We build the quad
> and hex QTrap rules via tensor products, but that's not how we order
> nodes.
> ---
> Roy
Thank you,
then how to solve the following problem.

I have an element that has n_nodes;
It has the same number of quadrature points that are at the nodal 
positions (trapezoidal rule is assumed).
How to find the node number that corresponds to a quadrature point?
Should I compare coordinates?

Thank you,
Michael.

Re: [Libmesh-users] question on Trapezoidal integration rule

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

On Wed, 7 Nov 2018, Michael Povolotskyi wrote:

> if I use the trapezoidal rule:
>
> QTrap qrule(dim)
>
> Can I assume the quadrature points are ordered in the same way at the element 
> nodes?

Doesn't look that way, in general, at first glance.  We build the quad
and hex QTrap rules via tensor products, but that's not how we order
nodes.
---
Roy

[Libmesh-users] question on Trapezoidal integration rule

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

Dear libmesh developers,

if I use the trapezoidal rule:

QTrap qrule(dim)

Can I assume the quadrature points are ordered in the same way at the 
element nodes?

Thank you,

Michael.

Re: [Libmesh-users] Quick question about libmesh: system.variable_type vs dof_map.variable_type?

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

Hi Roy,

Got it. Thank you so much!

Best,
Shawn

On Sun, Nov 4, 2018 at 5:58 AM Stogner, Roy H <roy...@ic...>
wrote:

>
> On Sat, 3 Nov 2018, Yuxiang Wang wrote:
>
> > Sorry again to spam with a novice question. As I read through libmesh
> > examples, I noticed that sometimes we get the variable_type from
> > dof_map.variable_type (where dof_map is a DofMap instance) and sometimes
> we
> > get from system.variable_type (where system is, for example, a
> > LinearImplicitSystem instance).
> >
> > Is there any difference between system.variable_type and
> > dof_map.variable_type? Are they pulling the same information?
>
> They're pulling the same information.
> ---
> Roy
>


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

Re: [Libmesh-users] Store the finite element matrices for post-processing of stresses/strains?

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

Hi David,

Thank you for your comment! Yes - for example, I found that ABAQUS mostly
use the average nodal values. In the mean time, it can have the user to
suppress averaging between domains of different materials or section
properties (e.g. cross-section of a beam/shell), in order to emphasize on
the discontinuities at the boundaries. I guess after all this is left to
user interpretation based on the actual physical problem being solved.

One of my friend also said that this is "not that critical" because it does
not impact the solving process. For the most important analyses, he usually
just uses quadrature stresses as the most "ground-truth" answer :)

Thank you again for the help David!

Best,
Shawn

On Sun, Nov 4, 2018 at 6:34 AM David Knezevic <dav...@ak...>
wrote:

> On Sat, Nov 3, 2018 at 11:46 PM Yuxiang Wang <yw...@vi...> wrote:
>
>> Hi David,
>>
>> Greetings again.
>>
>> I found the answer from this link
>> <http://www.joinville.udesc.br/portal/professores/pablo/materiais/CIL26_0053.pdf>
>> and yes, we will get a globally smooth distribution at the cost of much
>> more computation.
>>
>
> That looks like a good reference. I think the most common approach to
> getting a smooth stress plot is to compute element-based stress and then do
> nodal averaging. But there are obviously many different approaches that one
> can use.
>
> David
>
>
>
>>
>> On Fri, Nov 2, 2018 at 10:01 PM Yuxiang Wang <yw...@vi...> wrote:
>>
>>> Hi David,
>>>
>>> Thank you for your response! Yes I was indeed reading your example and
>>> thinking about it. Your insight helps a lot. Thank you so much!
>>>
>>> Since we are discussing on this topic, if you wouldn't mind, I'd like to
>>> ask one more question:
>>>
>>> During the L2 projection process, we solved the system per element. If
>>> we assemble the element mass matrices and the RHS to a global one with all
>>> elements, and similarly solve the global equation, would we be getting a
>>> globally smoothed stress distribution?
>>>
>>> I'd really love to hear your comment. Thanks in advance!
>>>
>>> Best,
>>> Shawn
>>>
>>> On Fri, Nov 2, 2018 at 7:27 PM David Knezevic <
>>> dav...@ak...> wrote:
>>>
>>>> On Fri, Nov 2, 2018 at 9:50 PM Yuxiang Wang <yw...@vi...> wrote:
>>>>
>>>>> Dear all,
>>>>>
>>>>> I have been using the matrix notation (following the Bathe textbook)
>>>>> with
>>>>> libmesh recently to code finite element solutions for basic continuum
>>>>> structural elements. I realized that I need to compute those matrices
>>>>> twice
>>>>> - once when I assemble the global stiffness matrix, and another time
>>>>> when I
>>>>> get the solutions and need to post-process to compute the
>>>>> stresses/strains.
>>>>> Being curious, do we have any best practices (or special
>>>>> considerations) to
>>>>> save those matrices so that we don't have to compute them twice? For
>>>>> example, should I just create a huge vector of DenseMatrix, and each
>>>>> matrix
>>>>> for each quadrature point? Or that libmesh has some tools for this book
>>>>> keeping?
>>>>>
>>>>> Explanation for what I meant by "matrix notation" and what are those
>>>>> finite
>>>>> element matrices: for example, for each quadrature point I have an
>>>>> interpolation matrix [H], a strain-displacement matrix [B] (constructed
>>>>> from derivatives of [H]). With the constitutive tensor matrix being
>>>>> [C], we
>>>>> can easily get the element stiffness matrix by [K] =
>>>>> integrate([B]^T[C][B]). After the solution [uhat] is obtained, we can
>>>>> also
>>>>> get the strain and stress by again [strain] = [B][uhat] or [stress] =
>>>>> [C][B][uhat].
>>>>>
>>>>> Just wondering whether anyone else has done this before and whether
>>>>> there
>>>>> would be a better practice to do it more elegantly. Thanks :)
>>>>>
>>>>> Best,
>>>>> Shawn
>>>>
>>>>
>>>> See systems_of_equations_ex6 for an example of computing stresses in
>>>> libMesh. In that case we do not store and re-use the per-element strain and
>>>> stress matrices (C and B in your notation) You could store and reuse that
>>>> data, but it's not clear to me whether it'd be worth it. You'd presumably
>>>> get some speedup in the evaluation of the stress at the cost of extra
>>>> storage, so you'd have to decide if you want that or not (in my experience
>>>> I would not want that since memory is usually more of a limiting factor
>>>> than speed).
>>>>
>>>> Best,
>>>> David
>>>>
>>>
>>>
>>> --
>>> Yuxiang "Shawn" Wang, PhD
>>> yw...@vi...
>>> +1 (434) 284-0836
>>>
>>
>>
>> --
>> Yuxiang "Shawn" Wang, PhD
>> yw...@vi...
>> +1 (434) 284-0836
>>
>

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

Re: [Libmesh-users] Store the finite element matrices for post-processing of stresses/strains?

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

On Sat, Nov 3, 2018 at 11:46 PM Yuxiang Wang <yw...@vi...> wrote:

> Hi David,
>
> Greetings again.
>
> I found the answer from this link
> <http://www.joinville.udesc.br/portal/professores/pablo/materiais/CIL26_0053.pdf>
> and yes, we will get a globally smooth distribution at the cost of much
> more computation.
>

That looks like a good reference. I think the most common approach to
getting a smooth stress plot is to compute element-based stress and then do
nodal averaging. But there are obviously many different approaches that one
can use.

David



>
> On Fri, Nov 2, 2018 at 10:01 PM Yuxiang Wang <yw...@vi...> wrote:
>
>> Hi David,
>>
>> Thank you for your response! Yes I was indeed reading your example and
>> thinking about it. Your insight helps a lot. Thank you so much!
>>
>> Since we are discussing on this topic, if you wouldn't mind, I'd like to
>> ask one more question:
>>
>> During the L2 projection process, we solved the system per element. If we
>> assemble the element mass matrices and the RHS to a global one with all
>> elements, and similarly solve the global equation, would we be getting a
>> globally smoothed stress distribution?
>>
>> I'd really love to hear your comment. Thanks in advance!
>>
>> Best,
>> Shawn
>>
>> On Fri, Nov 2, 2018 at 7:27 PM David Knezevic <dav...@ak...>
>> wrote:
>>
>>> On Fri, Nov 2, 2018 at 9:50 PM Yuxiang Wang <yw...@vi...> wrote:
>>>
>>>> Dear all,
>>>>
>>>> I have been using the matrix notation (following the Bathe textbook)
>>>> with
>>>> libmesh recently to code finite element solutions for basic continuum
>>>> structural elements. I realized that I need to compute those matrices
>>>> twice
>>>> - once when I assemble the global stiffness matrix, and another time
>>>> when I
>>>> get the solutions and need to post-process to compute the
>>>> stresses/strains.
>>>> Being curious, do we have any best practices (or special
>>>> considerations) to
>>>> save those matrices so that we don't have to compute them twice? For
>>>> example, should I just create a huge vector of DenseMatrix, and each
>>>> matrix
>>>> for each quadrature point? Or that libmesh has some tools for this book
>>>> keeping?
>>>>
>>>> Explanation for what I meant by "matrix notation" and what are those
>>>> finite
>>>> element matrices: for example, for each quadrature point I have an
>>>> interpolation matrix [H], a strain-displacement matrix [B] (constructed
>>>> from derivatives of [H]). With the constitutive tensor matrix being
>>>> [C], we
>>>> can easily get the element stiffness matrix by [K] =
>>>> integrate([B]^T[C][B]). After the solution [uhat] is obtained, we can
>>>> also
>>>> get the strain and stress by again [strain] = [B][uhat] or [stress] =
>>>> [C][B][uhat].
>>>>
>>>> Just wondering whether anyone else has done this before and whether
>>>> there
>>>> would be a better practice to do it more elegantly. Thanks :)
>>>>
>>>> Best,
>>>> Shawn
>>>
>>>
>>> See systems_of_equations_ex6 for an example of computing stresses in
>>> libMesh. In that case we do not store and re-use the per-element strain and
>>> stress matrices (C and B in your notation) You could store and reuse that
>>> data, but it's not clear to me whether it'd be worth it. You'd presumably
>>> get some speedup in the evaluation of the stress at the cost of extra
>>> storage, so you'd have to decide if you want that or not (in my experience
>>> I would not want that since memory is usually more of a limiting factor
>>> than speed).
>>>
>>> Best,
>>> David
>>>
>>
>>
>> --
>> Yuxiang "Shawn" Wang, PhD
>> yw...@vi...
>> +1 (434) 284-0836
>>
>
>
> --
> Yuxiang "Shawn" Wang, PhD
> yw...@vi...
> +1 (434) 284-0836
>