[deal.II] Re: Difficulty setting manifold with Opencascade using gmsh mesh + salome STEP (or IGES) - Step-54

[Bruno Blais]

Seems the image did not work,


On Wednesday, April 1, 2020 at 9:34:35 AM UTC-4, Bruno Blais wrote:
>
> So it seems part of my problem relate to the topology I was trying to 
> adapt OR to the staring complex mesh I was using.
> I have managed to make it work starting with a trivial mesh and using a 
> more simple IGES CAD. I will work on complexifying it from there. I will 
> try to keep the information here as it seems that feature has not been used 
> extensively :)!.
> What I have working now is this simple gif.
>
> [image: myimage.gif]
>
>
> On Tuesday, March 31, 2020 at 9:19:23 AM UTC-4, Bruno Blais wrote:
>>
>> Dear all,
>> I hope you are well in these uncertain times.
>>
>> I have been trying to use the OpenCASCADE library to set-up my manifolds. 
>> The goal of this is to be able to do CFD simulations in complex geometry, 
>> but use an initial very coarse mesh made with GMSH and use the dynamics 
>> mesh adaptation to adapt the mesh while having the adaptation be CAD-aware. 
>> The geometry is complex, so fixing the manifold analytically appears to be 
>> a bad idea.
>> I started from step-54, which compiles and works well on my machines.
>> However, I am unable to "make it work" with my own generated Salome STEP 
>> and gmsh MSH files.
>> My CAD is in metric, so I adjusted the scaling to use a scaling of 1. My 
>> mesh also contains numerous wires, so I disabled that aspect and only 
>> focused on refining with the faces which I label by looping over them.
>> However, no matter what occurs I get an error thrown:
>>
>>
>> 
>> An error occurred in line <114> of file 
>>  
>> in function
>> dealii::Point 
>> dealii::OpenCASCADE::NormalProjectionManifold> spacedim>::project_to_manifold(const dealii::ArrayView> dealii::Point >&, const dealii::Point&) const [with int 
>> dim = 2; int spacedim = 3]
>> The violated condition was: 
>> closest_point(sh, surrounding_points[i], tolerance) 
>> .distance(surrounding_points[i]) < std::max(tolerance * 
>> surrounding_points[i].norm(), tolerance)
>> Additional information: 
>> The point [ 0.1 0 0 ] is not on the manifold.
>>
>> Stacktrace:
>> ---
>> #0  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
>> dealii::OpenCASCADE::NormalProjectionManifold<2, 
>> 3>::project_to_manifold(dealii::ArrayView const, 
>> dealii::MemorySpace::Host> const&, dealii::Point<3, double> const&) const
>> #1  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
>> dealii::FlatManifold<2, 
>> 3>::get_new_point(dealii::ArrayView const, 
>> dealii::MemorySpace::Host> const&, dealii::ArrayView> dealii::MemorySpace::Host> const&) const
>> #2  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
>> dealii::Manifold<2, 
>> 3>::get_new_point_on_line(dealii::TriaIterator> 3> > const&) const
>> #3  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
>> #4  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
>> #5  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
>> dealii::TriaAccessor<1, 2, 3>::center(bool, bool) const
>> #6  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
>> dealii::Triangulation<2, 3>::DistortedCellList 
>> dealii::internal::TriangulationImplementation::Implementation::execute_refinement<3>(dealii::Triangulation<2,
>>  
>> 3>&, bool)
>> #7  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
>> dealii::Triangulation<2, 3>::execute_refinement()
>> #8  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
>> dealii::Triangulation<2, 3>::execute_coarsening_and_refinement()
>> #9  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
>> dealii::Triangulation<2, 3>::refine_global(unsigned int)
>> #10  ./step-54: Step54::TriangulationOnCAD::refine_mesh()
>> #11  ./step-54: Step54::TriangulationOnCAD::run()
>> #12  ./step-54: main
>> 
>>
>> I have tried some elements :
>> - Playing with the tolerance of the CAD
>> - Writing my CAD from the OpenCASCADE STEP object and then opening it 
>> again in Salome to see that it is interpreted correctly +
>>
>> However, I always seem to face an issue where my CAD does not appear to 
>> generate a valid manifold. I have joined my code as well as the salome HDF 
>> file and gmsh GEO file used to generate the CAD and meshes respectively.
>> I must admit I am quite at lost here as to why both don't coincide. Any 
>> advices on where to proceed from there?
>> Best :)!
>> Bruno
>>
>>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To 

[deal.II] Re: Difficulty setting manifold with Opencascade using gmsh mesh + salome STEP (or IGES) - Step-54

[Bruno Blais]

So it seems part of my problem relate to the topology I was trying to adapt 
OR to the staring complex mesh I was using.
I have managed to make it work starting with a trivial mesh and using a 
more simple IGES CAD. I will work on complexifying it from there. I will 
try to keep the information here as it seems that feature has not been used 
extensively :)!.
What I have working now is this simple gif.

[image: myimage.gif]


On Tuesday, March 31, 2020 at 9:19:23 AM UTC-4, Bruno Blais wrote:
>
> Dear all,
> I hope you are well in these uncertain times.
>
> I have been trying to use the OpenCASCADE library to set-up my manifolds. 
> The goal of this is to be able to do CFD simulations in complex geometry, 
> but use an initial very coarse mesh made with GMSH and use the dynamics 
> mesh adaptation to adapt the mesh while having the adaptation be CAD-aware. 
> The geometry is complex, so fixing the manifold analytically appears to be 
> a bad idea.
> I started from step-54, which compiles and works well on my machines.
> However, I am unable to "make it work" with my own generated Salome STEP 
> and gmsh MSH files.
> My CAD is in metric, so I adjusted the scaling to use a scaling of 1. My 
> mesh also contains numerous wires, so I disabled that aspect and only 
> focused on refining with the faces which I label by looping over them.
> However, no matter what occurs I get an error thrown:
>
>
> 
> An error occurred in line <114> of file 
>  
> in function
> dealii::Point 
> dealii::OpenCASCADE::NormalProjectionManifold spacedim>::project_to_manifold(const dealii::ArrayView dealii::Point >&, const dealii::Point&) const [with int 
> dim = 2; int spacedim = 3]
> The violated condition was: 
> closest_point(sh, surrounding_points[i], tolerance) 
> .distance(surrounding_points[i]) < std::max(tolerance * 
> surrounding_points[i].norm(), tolerance)
> Additional information: 
> The point [ 0.1 0 0 ] is not on the manifold.
>
> Stacktrace:
> ---
> #0  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
> dealii::OpenCASCADE::NormalProjectionManifold<2, 
> 3>::project_to_manifold(dealii::ArrayView const, 
> dealii::MemorySpace::Host> const&, dealii::Point<3, double> const&) const
> #1  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
> dealii::FlatManifold<2, 
> 3>::get_new_point(dealii::ArrayView const, 
> dealii::MemorySpace::Host> const&, dealii::ArrayView dealii::MemorySpace::Host> const&) const
> #2  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
> dealii::Manifold<2, 
> 3>::get_new_point_on_line(dealii::TriaIterator 3> > const&) const
> #3  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
> #4  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
> #5  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
> dealii::TriaAccessor<1, 2, 3>::center(bool, bool) const
> #6  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
> dealii::Triangulation<2, 3>::DistortedCellList 
> dealii::internal::TriangulationImplementation::Implementation::execute_refinement<3>(dealii::Triangulation<2,
>  
> 3>&, bool)
> #7  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
> dealii::Triangulation<2, 3>::execute_refinement()
> #8  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
> dealii::Triangulation<2, 3>::execute_coarsening_and_refinement()
> #9  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
> dealii::Triangulation<2, 3>::refine_global(unsigned int)
> #10  ./step-54: Step54::TriangulationOnCAD::refine_mesh()
> #11  ./step-54: Step54::TriangulationOnCAD::run()
> #12  ./step-54: main
> 
>
> I have tried some elements :
> - Playing with the tolerance of the CAD
> - Writing my CAD from the OpenCASCADE STEP object and then opening it 
> again in Salome to see that it is interpreted correctly +
>
> However, I always seem to face an issue where my CAD does not appear to 
> generate a valid manifold. I have joined my code as well as the salome HDF 
> file and gmsh GEO file used to generate the CAD and meshes respectively.
> I must admit I am quite at lost here as to why both don't coincide. Any 
> advices on where to proceed from there?
> Best :)!
> Bruno
>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/a239784b-954a-4860-ab3a-acac05e4c709%40googlegroups.com.

Re: [deal.II] step-42 now available

[Muhammad Mashhood]

Hi Mac,
 I really appreciate your concern and thank you for nice and 
brief summary regarding classical plasticity algorithm. 
I would like to highlight here that in Prisms-plasticity this algorithm is 
also implemented but I checked with its continuum plasticity application, 
which so far is not updated and hence also the code is not working for 
continuum plasticity part. As per communication from developers they will 
update and make it run in next update.
Secondly I am preferably using the step-42. Personally for me it seems bit 
simpler and I already coupled it with other physics for multi-physics 
problem.


On Wednesday, April 1, 2020 at 11:26:44 AM UTC+2, mac wrote:
>
> Hi all 
>
> Classical plasticity is underpinned by return mapping algorithms that 
> operate at the level of the quadrature point. One has a global predictor 
> for the displacement field where we assume frozen plastic flow. Then at 
> each quadrature point you compute a trial strain, and hence a trial stress. 
> The trial stress is used to test for yield at the quadrature point. You 
> then compute the plastic multiplier, and hence the plastic strain, to 
> return to the yield surface. This is an iterative process at the global and 
> possibly local level.
>
No doubt, you are very right that the plastic multiplier is evaluated in 
iterative manner and this very plastic multiplier is then a key to plastic 
strain evaluation. But here I would like to highlight that in the step-42 
seemingly a bit different method is used to solve elastoplastic problem. As 
far as I could understand the algorithm in hand, here after evaluating the 
trial stress in terms of Von-Mises equivalent stress and also when it is 
higher than the yield stress, the stiffness tensor is evaluated as function 
of current trial stress. After that this very stiffness tensor is used to 
evaluate displacement and strain increments and convergence of Newton 
algorithm is checked. So in this way the system is solved iteratively and 
after convergence the displacement and the total strain is present as 
output. 
As the total strain is the output that is why I need here some way to 
evaluate the plastic and elastic strain (any one of them so that other can 
be evaluated from additive decomposition). If you want, you can check this 
interesting and simple algorithm in a deal.ii step-42 code or this nice 
article of "Goal oriented error estimation and mesh adaptivity in 3d 
elastoplasticity problems" doi: 
https://doi-org.proxy.bnl.lu/10.1007/s10704-016-0113-y 

I hope I explained it well. Thanks again for showing interest and your 
concern. 

>
> All of this can and has been done in deal.ii. Wolfgang has suggested using 
> CellDataStorage which will allow you to handle the quadrature point data. 
> For implementations have a look at 
> https://github.com/prisms-center/plasticity . 
>
> You might also want to have a look at some of the extensive literature 
> e.g. https://onlinelibrary.wiley.com/doi/book/10.1002/9780470694626 Here 
> they explain the algorithm in detail.
>
> Best
> A
>
> On 30 Mar 2020, at 19:09, Muhammad Mashhood  > wrote:
>
>
>
> On Monday, February 24, 2020 at 12:42:17 AM UTC+1, Wolfgang Bangerth wrote:
>>
>>
>> > In my question I meant whether is it possible to evaluated plastic 
>> strain 
>> > component from currently implemented plasticity algorithm as a further 
>> > development of this code? 
>>
>> I am pretty sure that it is, by just computing the difference between the 
>> elastic stress (C eps(u)) and the actual stress computed. In fact, for 
>> the 
>> current situation, the actual stress computed equal to the elastic stress 
>> where it is less than the yield stress (and so the plastic strain is 
>> zero), 
>> and it is simply a fraction of the elastic stress where it exceeds the 
>> yield 
>> stress. Once you have the plastic stress, you can compute the plastic 
>> strain 
>> by multiplying it by C^{-1}. 
>> Thank you very much Prof. Wolfgang for your suggestion and I hope you are 
>> doing all well. As far as I understood now from your suggestion is that on 
>> the base of additive decomposition, I can subtract elastic stress (C 
>> eps(u)) from the actual stress computed to have the plastic stress part. 
>> But here I just noticed in the program that we also need (eps(u)) or the 
>> elastic strain part to have elastic stress part's value when the domain is 
>> in the plastic region.
>
> In one dimension problem there is I think no problem because from 
> stress strain curve one knows already the elastic stress value (which is 
> yield stress) but here we have multidimensional problem i.e. stress and 
> strain as a 2nd order tensor. Therefore I am thinking now that how this 
> elastic strain tensor can be evaluated to have elastic stress part to be 
> subtracted from total calculated stress. 
> One of the idea I have so far is to use the calculated "strain_tensor" to 
> find the elastic stress as  "elastic_stress_tensor = 
> 

Re: [deal.II] step-42 now available

[Andrew McBride]

Hi all 

Classical plasticity is underpinned by return mapping algorithms that operate 
at the level of the quadrature point. One has a global predictor for the 
displacement field where we assume frozen plastic flow. Then at each quadrature 
point you compute a trial strain, and hence a trial stress. The trial stress is 
used to test for yield at the quadrature point. You then compute the plastic 
multiplier, and hence the plastic strain, to return to the yield surface. This 
is an iterative process at the global and possibly local level.

All of this can and has been done in deal.ii. Wolfgang has suggested using 
CellDataStorage which will allow you to handle the quadrature point data. For 
implementations have a look at https://github.com/prisms-center/plasticity 
 . 

You might also want to have a look at some of the extensive literature e.g. 
https://onlinelibrary.wiley.com/doi/book/10.1002/9780470694626 
 Here they 
explain the algorithm in detail.

Best
A

> On 30 Mar 2020, at 19:09, Muhammad Mashhood  wrote:
> 
> 
> 
> On Monday, February 24, 2020 at 12:42:17 AM UTC+1, Wolfgang Bangerth wrote:
> 
> > In my question I meant whether is it possible to evaluated plastic strain 
> > component from currently implemented plasticity algorithm as a further 
> > development of this code? 
> 
> I am pretty sure that it is, by just computing the difference between the 
> elastic stress (C eps(u)) and the actual stress computed. In fact, for the 
> current situation, the actual stress computed equal to the elastic stress 
> where it is less than the yield stress (and so the plastic strain is zero), 
> and it is simply a fraction of the elastic stress where it exceeds the yield 
> stress. Once you have the plastic stress, you can compute the plastic strain 
> by multiplying it by C^{-1}. 
> Thank you very much Prof. Wolfgang for your suggestion and I hope you are 
> doing all well. As far as I understood now from your suggestion is that on 
> the base of additive decomposition, I can subtract elastic stress (C eps(u)) 
> from the actual stress computed to have the plastic stress part. But here I 
> just noticed in the program that we also need (eps(u)) or the elastic strain 
> part to have elastic stress part's value when the domain is in the plastic 
> region.
> In one dimension problem there is I think no problem because from stress 
> strain curve one knows already the elastic stress value (which is yield 
> stress) but here we have multidimensional problem i.e. stress and strain as a 
> 2nd order tensor. Therefore I am thinking now that how this elastic strain 
> tensor can be evaluated to have elastic stress part to be subtracted from 
> total calculated stress. 
> One of the idea I have so far is to use the calculated "strain_tensor" to 
> find the elastic stress as  "elastic_stress_tensor = 
> (stress_strain_tensor_kappa + stress_strain_tensor_mu) * strain_tensor". But 
> again the point is that the "strain_tensor" here is the total strain rather 
> than elastic strain component.  
> 
> Kindly correct me if I misunderstood your suggestion at any point or if there 
> is an alternate approach possible. Thank you again in advance. 
> Stay healthy stay safe!
> 
> 
> > Then as a step further I would be trying to store this plastic strain in 
> > cells 
> > or Gauss points along with the modified yield strength (due to isotropic 
> > hardening) so that history of loading is stored too in the domain. 
> 
> I imagine that that, too, can be done. I'm not an expert in plasticity, but I 
> see no fundamental reasons why what you want to do should not be possible. 
> There are also classes CellDataStorage and TransferableQuadraturePointData 
> and 
> parallel::distributed::ContinuousQuadratureDataTransfer that can help you 
> with 
> storing information at quadrature points. 
> 
> Best 
>   W. 
> 
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>  
> 
> 
> -- 
> The deal.II project is located at http://www.dealii.org/ 
> 
> For mailing list/forum options, see 
> https://groups.google.com/d/forum/dealii?hl=en 
> 
> --- 
> You received this message because you are subscribed to the Google Groups 
> "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an 
> email to dealii+unsubscr...@googlegroups.com 
> .
> To view this discussion on the web visit 
> https://groups.google.com/d/msgid/dealii/9a59d3e4-e11d-4f6e-af92-c34c85fe9b57%40googlegroups.com
>  
>