RE: [deal.II] How to get the shape_value at the quadrature point in the reference cell?

[Michael Li]

Dr. Bangerth, Thank you for your hints. I tried using shape_value, but it seems  finite_element_output is a protected member. The error shows as follows: error: ‘dealii::internal::FEValuesImplementation::FiniteElementRelatedData<2, 2> dealii::FEValuesBase<2, 2>::finite_element_output’ is protected within this context   59 |   double shape_value = fe_values.finite_element_output.shape_value(i, q_index);   ^error: ‘class dealii::internal::FEValuesImplementation::FiniteElementRelatedData<2, 2>’ has no member named ‘shape_value’; did you mean ‘shape_values’?  By the way, I also want to get the shape_grad of the reference cell to compare with some other code implementation regarding shape function.  Based on Thong’s method, the following codes gives me the desired results for a Q1 element:  const FE_Q<2> Q_2d_element(1);  const Point<2> center(1, 1);  for(unsigned int i=0; i  {std::cout << " shape value = " << Q_2d_element.shape_value(i, center);std::cout << " shape grad = "  << Q_2d_element.shape_grad(i, center) << std::endl;  } I wonder if there are other ways to access the shape_grad on the ref cell directly. Do I need to use fe_values.inverse_jacobian(q_index) to convert from shape_grad on the real cell to the reference cell?  Best,Michael From: Wolfgang BangerthSent: Sunday, September 12, 2021 5:13 PMTo: dealii@googlegroups.comSubject: Re: [deal.II] How to get the shape_value at the quadrature point in the reference cell? On 9/10/21 4:48 PM, Michael Lee wrote:> > I want to require the value of a shape function at some quadrature point. But > the following code throws an error saying finite_element_output is protected.> double shape_value = fe_values.finite_element_output.shape_values(i, q_index);>  You already solved the problem, but just for completeness: You get the error because you try to access a 'protected' member variable called 'shape_values' (plural). What you want is to use the 'public' member function 'shape_value' (singular). Best  W. -- Wolfgang Bangerth  email: bange...@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/8cd6a2e7-5ad5-fbf5-3240-3237352335a9%40colostate.edu. 



-- 
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/8ACD76BA-4163-4FFB-832E-01D37DCF18C0%40hxcore.ol.

RE: [deal.II] Re: How to get the shape_value at the quadrature point in the reference cell?

[Michael Li]

Hi Dear Thang,Thank you for you code snippet! It is simple and works well and is what I want. The link helps me understand the shape function in deal.ii. Have a good weekend!Michael  From: Thang W PhamSent: Saturday, September 11, 2021 12:20 AMTo: deal.II User GroupSubject: [deal.II] Re: How to get the shape_value at the quadrature point in the reference cell? Hi Michael,I think you should have a look at FE_Q (Lagrange elements) at this link. In deal.II, the reference cell's domain is [0,1]^2 for 2d element. So, the point you are about to evaluate (to get the value of phi_i equal 0.25)  is (0.5, 0.5).You can refer to this snippet:/*  const FE_Q<2> Q_2d_element(1);  const Point<2> center(0.5, 0.5);  for(int i=0; istd::cout << Q_2d_element.shape_value(i, center) << std::endl;  }*///output 0.25 0.25 0.25 0.25Also  have a look at step-3 to understand how FE_Q works in real problem since the shape functions will be evaluated at real cell coordinates. Good luck. Regards,T.P Vào lúc 07:48:35 UTC+9 ngày Thứ Bảy, 11 tháng 9, 2021, lian...@gmail.com đã viết:Deal deal.ii comunity, I want to require the value of a shape function at some quadrature point. But the following code throws an error saying finite_element_output is protected.     double shape_value = fe_values.finite_element_output.shape_values(i, q_index); For a 2D bilinear element, I expect the shape values  at the cell center(kesi=0, eta=0) are 0.25 for all the 4 shape functions.              Is there other way to acquire this information? Best,Michael-- 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/b939b310-fd2b-46ee-83fa-59b46e6f9130n%40googlegroups.com. 



-- 
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/11C955E7-916F-420A-94AD-9479516E24F8%40hxcore.ol.

RE: [deal.II] How to watch the variable of BlockVector type?

[Michael Li]

Yeah! .block(0) is essential. Without this, the data obtained is garbage! Thanks for this great hint! From: Wolfgang BangerthSent: Wednesday, August 18, 2021 8:37 PMTo: dealii@googlegroups.comSubject: Re: [deal.II] How to watch the variable of BlockVector type? On 8/18/21 7:18 PM, Michael Lee wrote:> >            BlockVector & system_rhs;>            system_rhs.reinit(dofs_per_block);> > For normal vectors, I can use *((double(*)[10]) v.begin()), but it does not > work for BlockVector. I'm using VS code. Try *( (double(*)[10]) v.block(0).begin()) Best  W.  -- Wolfgang Bangerth  email: bange...@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/46c327fa-1d8c-0eb9-d03e-40e5c3611447%40colostate.edu. 



-- 
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/B0E16C6C-D172-4D48-BD6D-30EEC923A16A%40hxcore.ol.

RE: [deal.II] How to watch the variable of BlockVector type?

[Michael Li]

Sorry, it seems to work using *((double(*)[20])residual.begin()).   From: Michael LeeSent: Wednesday, August 18, 2021 7:18 PMTo: deal.II User GroupSubject: [deal.II] How to watch the variable of BlockVector type? I want to check if I have system_rhs assembled correctly. Is there an easy way to watch the contents of this BlockVector during debugging?           BlockVector & system_rhs;          system_rhs.reinit(dofs_per_block); For normal vectors, I can use *((double(*)[10]) v.begin()), but it does not work for BlockVector. I'm using VS code. Best,Michael-- 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/9db1415e-812b-4825-8f73-eb7ed3866c17n%40googlegroups.com. 



-- 
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/B3897C92-3FC3-4B18-88D2-185E7F2AD4D9%40hxcore.ol.

RE: [deal.II] Computing the solution gradient at the quadrature point on a face

[Michael Li]

Hi Jean-Pau, You’re always helpful. The deformation dependent loads are discussed in Wriggers, Wood and some other books and papers. I’m studying that. It seems the formulation is complicated. I believe that the auto-differentiation tools is a good approach to this nonlinear problem and I’ll dig into it later. For now I hope I can have a better understanding this topic and implement it in deal.II by hand then I will be comfortable with AD. Thank you!Michael  From: Jean-Paul PelteretSent: Sunday, August 8, 2021 12:16 PMTo: dealii@googlegroups.comSubject: Re: [deal.II] Computing the solution gradient at the quadrature point on a face Hi Michael, I’m glad that you’re on the right track now.  If you could give some advice or refer to some references on this topic, I would appreciate that very mush. Unfortunately, this is where pretty much I hit the limit of what I can say on this topic as it would conflict with the expectations of my employer. There’s a little by-line on one approach in this paperhttps://www.sciencedirect.com/science/article/abs/pii/S0997753814001508But as for other references, I’m hoping that someone else can give you some guidance there. (There’s at least one fairly common book in nonlinear mechanics that discusses the "linearisation of external virtual work”. Maybe a search in this direction might get you somewhere.) Another option would be to try the auto-differentiation tools to do all of the heavy lifting for you. With few restrictions, it would be compatible with any general traction load. At the very least and AD implementation it could help you validate any linearisation that you do by hand. You can take a look at steps 70 and 71 to get some idea as to what it might be able to do for you. Using the approach in step-71 is probably easiest to start off in your case. You could even limit it to just linearising the contribution that comes from this boundary / traction contribution. Best,Jean-Paul On 4. Aug 2021, at 23:22, Michael Li <lianxi...@gmail.com> wrote: Hi Jean-Paul, Yes, I set up the boundary check when taking care of the traction. It worked as expected using Physics::Transformations::nansons_formula(face_normal, F). To make a quick test, I did not linearize the load stiffness but just used the load at last Newton-Raphson step, so it is like a mixture of N-R and Picard iteration. The technique is not correct and there is some discrepancy compared with the reference. But the result is better than not considering the face area change at all. So in some sense, it tells me that the surface transform is correct. I’m working on how to linearize an arbitrary surface traction. I have found some references dealing with the normal(pressure) load but not arbitrary surface load (the force is not normal to the face in my case). If you could give some advice or refer to some references on this topic, I would appreciate that very mush. You’re totally right, it converges only if a very small load step is used. The residual grows up fast and the result is meaningless if the load increment is large. Thanks for your great help! Michael   From: Jean-Paul PelteretSent: Wednesday, August 4, 2021 12:19 AMTo: dealii@googlegroups.comSubject: Re: [deal.II] Computing the solution gradient at the quadrature point on a face Hi Michael, The one thing that stands out about the snippet of code that you shared is that there is no check to if (cell->face(face)->at_boundary() == true && )...Do you still have something like that? Otherwise, after a cursory glance, everything appears to be alright with this. It mimics what you’d do to get the deformation gradient at a QP in a cell.  To debug further, what happens when you apply a very small load? Maybe its not scaled appropriately for the stiffness of the material. Also, if you’ve now added a deformation dependent part to the load (i.e. the area transformation) then you might need to consistently linearise this in order to attain convergence for large loads. Best,Jean-PaulOn 2. Aug 2021, at 18:03, Michael Li <lianxi...@gmail.com> wrote: Hi Jean-Paul, Thank you for your hints.  I initialized a local variable solution_grads_u_face_total with fe_face_values_ref[u_fe].get_function_gradients to get the gradient of displacement at the quadrature point on the face (see codes below). But the gradient acquired is not sensible and the resultant deformation is totally wrong. void assemble_neumann_contribution_one_cell(…){    ….  const FEValuesExtractors::Vector _fe = data.solid->u_fe; std::vector2, dim,NumberType> > solution_grads_u_face_total(data.solid->qf_face.size()); scratch.fe_face_values_ref.reinit(cell, face); scratch.fe_face_values_ref[u_fe].get_function_gradients(scratch.solution_total,solution_grads_u_face_total);   for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face)   {   for (unsigned int f_q_po

RE: [deal.II] Debug multithread code

[Michael Li]

I put MultithreadInfo::set_thread_limit(1) in main function, but still cannot do step debug in the function assemble_neumann_contribution_one_cell. Am I doing something wrong?    From: Timo HeisterSent: Thursday, August 5, 2021 1:19 PMTo: dealii@googlegroups.comSubject: Re: [deal.II] Debug multithread code Take a look at set_thread_limit:https://www.dealii.org/current/doxygen/deal.II/classMultithreadInfo.htmlAs you can see, you can also control it using an environment variable, so you don't have to recompile. On Thu, Aug 5, 2021, 13:36 Michael Lee  wrote:Hi, I find it's difficult to debug multithread code on a multiple-core computer. Sometimes it does not step into a function (I'm using vscode + wsl). Does deal.II use all the cores available? For the sake of debug, how can I limit the cores used? Is there a way to debug the code just like series one? Regards,Michael-- 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/b123aaf2-8d32-4b43-b818-8d2ae9ffac77n%40googlegroups.com.-- 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/CAMRj59HvBtTTF1D3SnGHgbnub4nNN_6TC0aYPqHJRK3NYoqdVA%40mail.gmail.com. 



-- 
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/B6A6B8CA-9E60-4912-BC84-5706D0E86574%40hxcore.ol.

RE: [deal.II] Computing the solution gradient at the quadrature point on a face

[Michael Li]

Hi Jean-Paul, Yes, I set up the boundary check when taking care of the traction. It worked as expected using Physics::Transformations::nansons_formula(face_normal, F). To make a quick test, I did not linearize the load stiffness but just used the load at last Newton-Raphson step, so it is like a mixture of N-R and Picard iteration. The technique is not correct and there is some discrepancy compared with the reference. But the result is better than not considering the face area change at all. So in some sense, it tells me that the surface transform is correct. I’m working on how to linearize an arbitrary surface traction. I have found some references dealing with the normal(pressure) load but not arbitrary surface load (the force is not normal to the face in my case). If you could give some advice or refer to some references on this topic, I would appreciate that very mush. You’re totally right, it converges only if a very small load step is used. The residual grows up fast and the result is meaningless if the load increment is large. Thanks for your great help! Michael    From: Jean-Paul PelteretSent: Wednesday, August 4, 2021 12:19 AMTo: dealii@googlegroups.comSubject: Re: [deal.II] Computing the solution gradient at the quadrature point on a face Hi Michael, The one thing that stands out about the snippet of code that you shared is that there is no check to if (cell->face(face)->at_boundary() == true && )...Do you still have something like that? Otherwise, after a cursory glance, everything appears to be alright with this. It mimics what you’d do to get the deformation gradient at a QP in a cell.  To debug further, what happens when you apply a very small load? Maybe its not scaled appropriately for the stiffness of the material. Also, if you’ve now added a deformation dependent part to the load (i.e. the area transformation) then you might need to consistently linearise this in order to attain convergence for large loads. Best,Jean-PaulOn 2. Aug 2021, at 18:03, Michael Li <lianxi...@gmail.com> wrote: Hi Jean-Paul, Thank you for your hints.  I initialized a local variable solution_grads_u_face_total with fe_face_values_ref[u_fe].get_function_gradients to get the gradient of displacement at the quadrature point on the face (see codes below). But the gradient acquired is not sensible and the resultant deformation is totally wrong. void assemble_neumann_contribution_one_cell(…){    ….  const FEValuesExtractors::Vector _fe = data.solid->u_fe; std::vector2, dim,NumberType> > solution_grads_u_face_total(data.solid->qf_face.size()); scratch.fe_face_values_ref.reinit(cell, face); scratch.fe_face_values_ref[u_fe].get_function_gradients(scratch.solution_total,solution_grads_u_face_total);   for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face)   {   for (unsigned int f_q_point = 0; f_q_point < n_q_points_f; ++f_q_point)   {   const Tensor<2,dim,NumberType> _u = solution_grads_u_face_total[f_q_point];   const Tensor<2,dim,NumberType> F = Physics::Elasticity::Kinematics::F(grad_u);   .       }   ………..}………..} Can you tell me what’s wrong with the code above? Thank you!Michael  From: Jean-Paul PelteretSent: Sunday, August 1, 2021 1:46 PMTo: dealii@googlegroups.comSubject: Re: [deal.II] Computing the solution gradient at the quadrature point on a face Hi Michael, The problem here is that “scratch.solution_grads_u_total” is initialised with the cell FEValues, and so has n_cell_quadrature_points entries. You need to pass "scratch.fe_face_values_ref[u_fe].get_function_gradients(scratch.solution_total,  scratch.solution_grads_u_total);” something that has size n_face_quadrature_points. That’s what this error is informing you.As a side note, Physics::Transformations::nansons_formula() can compute the transformation of the normal vector + area change for you, should you like the traction direction to follow the moving boundary surface as well. Best,Jean-PaulOn 1. Aug 2021, at 20:05, Michael Lee <lianxi...@gmail.com> wrote: fe_face_values.get_function_gradients returns a rank 2 tensor as the solution is the displacement vector. Let me be more specific. In the code gallery " Quasi-Static Finite-Strain Compressible Elasticity (The deal.II Library: The 'Quasi-Static Finite-Strain Compressible Elasticity' code gallery program (dealii.org) ). It does not consider the follower force as described in the documentation:  "                 // We specify the traction in reference configuration.// For this problem, a defined total vertical force is applied// in the reference configuration.// The direction of the applied traction is assumed not to// evolve with the deformation of the domain." So I made a modification to consider the effect of the deformed area

RE: [deal.II] Computing the solution gradient at the quadrature point on a face

[Michael Li]

Hi Jean-Paul, Thank you for your hints.  I initialized a local variable solution_grads_u_face_total with fe_face_values_ref[u_fe].get_function_gradients to get the gradient of displacement at the quadrature point on the face (see codes below). But the gradient acquired is not sensible and the resultant deformation is totally wrong. void assemble_neumann_contribution_one_cell(…){    ….  const FEValuesExtractors::Vector _fe = data.solid->u_fe; std::vector2, dim,NumberType> > solution_grads_u_face_total(data.solid->qf_face.size()); scratch.fe_face_values_ref.reinit(cell, face); scratch.fe_face_values_ref[u_fe].get_function_gradients(scratch.solution_total,solution_grads_u_face_total);   for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face)   {   for (unsigned int f_q_point = 0; f_q_point < n_q_points_f; ++f_q_point)   {   const Tensor<2,dim,NumberType> _u = solution_grads_u_face_total[f_q_point];   const Tensor<2,dim,NumberType> F = Physics::Elasticity::Kinematics::F(grad_u);   .       }   ………..}………..} Can you tell me what’s wrong with the code above? Thank you!Michael   From: Jean-Paul PelteretSent: Sunday, August 1, 2021 1:46 PMTo: dealii@googlegroups.comSubject: Re: [deal.II] Computing the solution gradient at the quadrature point on a face Hi Michael, The problem here is that “scratch.solution_grads_u_total” is initialised with the cell FEValues, and so has n_cell_quadrature_points entries. You need to pass "scratch.fe_face_values_ref[u_fe].get_function_gradients(scratch.solution_total,  scratch.solution_grads_u_total);” something that has size n_face_quadrature_points. That’s what this error is informing you.As a side note, Physics::Transformations::nansons_formula() can compute the transformation of the normal vector + area change for you, should you like the traction direction to follow the moving boundary surface as well. Best,Jean-PaulOn 1. Aug 2021, at 20:05, Michael Lee  wrote: fe_face_values.get_function_gradients returns a rank 2 tensor as the solution is the displacement vector. Let me be more specific. In the code gallery " Quasi-Static Finite-Strain Compressible Elasticity (The deal.II Library: The 'Quasi-Static Finite-Strain Compressible Elasticity' code gallery program (dealii.org) ). It does not consider the follower force as described in the documentation:  "                 // We specify the traction in reference configuration.// For this problem, a defined total vertical force is applied// in the reference configuration.// The direction of the applied traction is assumed not to// evolve with the deformation of the domain." So I made a modification to consider the effect of the deformed area on the traction as follows according to the formula da/dA = J sqrt(N C^-1 N) .     voidassemble_neumann_contribution_one_cell(const typename DoFHandler::active_cell_iterator ,   ScratchData_ASM ,   PerTaskData_ASM ){  // Aliases for data referenced from the Solid class  const unsigned int _q_points_f = data.solid->n_q_points_f;  const unsigned int _per_cell = data.solid->dofs_per_cell;  const Time  = data.solid->time;  const FESystem  = data.solid->fe;  const unsigned int _dof = data.solid->u_dof;  const FEValuesExtractors::Vector _fe = data.solid->u_fe;       // Next we assemble the Neumann contribution. We first check to see it the  // cell face exists on a boundary on which a traction is applied and add  // the contribution if this is the case.  for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face)if (cell->face(face)->at_boundary() == true&& cell->face(face)->boundary_id() == 11)  {scratch.fe_face_values_ref.reinit(cell, face);scratch.fe_face_values_ref[u_fe].get_function_gradients(scratch.solution_total,scratch.solution_grads_u_total); for (unsigned int f_q_point = 0; f_q_point < n_q_points_f; ++f_q_point)  {                // Note that the contributions to the right hand side vector we// compute here only exist in the displacement components of the// vector.const double time_ramp = (time.current() / time.end());const double magnitude  = 0.5*1.95*1000*0.32*0.32*time_ramp; // (Total force) / (RHS surface area)Tensor<1,dim> dir;dir[0] = 1.0;Tensor<1,dim> face_normal = scratch.fe_face_values_ref.normal_vector(f_q_point); const Tensor<2,dim,NumberType> _u = scratch.solution_grads_u_total[f_q_point];const Tensor<2,dim,NumberType> F = 

RE: [deal.II] Matrix Multiplication

[Michael Li]

Hi Jean-Paul,Thanks for confirming that. I have made the Kirchhoff model with your suggestion. The Physics::Transformations is so convenient for elasticity problems.  In case someone wonders how I did it, I enclose the codes here. SymmetricTensor<4,dim,NumberType>get_Jc(const Tensor< 2, dim, NumberType >  &F) const{  SymmetricTensor<4,dim, NumberType> D = get_stress_strain_tensor(lambda, mu);   return Physics::Transformations::Contravariant::push_forward(D, F);} SymmetricTensor<2,dim,NumberType>get_tau(const Tensor< 2, dim, NumberType >  &F){  // Second Piola Kirchhoff stress  SymmetricTensor<2,dim,NumberType> S = get_stress_strain_tensor(lambda, mu) *Physics::Elasticity::Kinematics::E(F);   return Physics::Transformations::Contravariant::push_forward(S, F);} I just undefine the variable to circumvent the Sacado error.  #define ENABLE_SACADO_FORMULATION   Thanks,Michael   From: Jean-Paul PelteretSent: Monday, July 26, 2021 12:46 AMTo: dealii@googlegroups.comSubject: Re: [deal.II] Matrix Multiplication Hi Michael, I revisited the code gallery example “The deal.II Library: The 'Quasi-Static Finite-Strain Compressible Elasticity' code gallery program (dealii.org)”. It seems I just need to modify two functions, i.e., getJc(det_F, b_bar) and get_tau(det_F, b_bar) if I want to implement the St. Vienant Kirchhoff material. Is that correct? Yes, that’s correct. You need only supply a new definition of the Kirchhoff stress and its linearisation in order to implement a new constitutive law. The child namespaces of Physics::Transformations can be of some assistance if you choose a different parameterisation (i.e. use a different strain measure leading to different natural stress measure) and need to push-forward the result into the current configuration. You also don’t need to adhere to the volumetric/isochoric split that is used in the code gallery example.  With regard to the error message that you’re seeing, this was a bug that was fixed in the upcoming 9.3 release.https://github.com/dealii/dealii/pull/12217The warning messages also relate to some functions that are going to be removed later, so we need to update the program to use the function which supersedes it. Thanks for the letting us know about it! Best,Jean-Paul On 22. Jul 2021, at 23:54, Michael Li <lianxi...@gmail.com> wrote: Hi Jean-Paul, I revisited the code gallery example “The deal.II Library: The 'Quasi-Static Finite-Strain Compressible Elasticity' code gallery program (dealii.org)”. It seems I just need to modify two functions, i.e., getJc(det_F, b_bar) and get_tau(det_F, b_bar) if I want to implement the St. Vienant Kirchhoff material. Is that correct? I compiled cook_membrane.cc a couple of days ago, but for no reason, it fails giving the following errors: /home/michael/dealii-9.2.0/examples/code-gallery/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc:1815:5:   required from here/home/michael/share/trilinos/include/Sacado_Fad_Exp_Ops.hpp:775:1: error: no type named ‘type’ in ‘struct Sacado::mpl::disable_if, Sacado::Fad::Exp::MultiplicationOp >, double, false, true, Sacado::Fad::Exp::ExprSpecDefault, Sacado::Fad::Exp::PowerImpl::Scalar>, double, false, true, Sacado::Fad::Exp::ExprSpecDefault> >’ /usr/local/include/deal.II/physics/elasticity/kinematics.h:294:47: error: no match for ‘operator*’ (operand types are ‘Sacado::Fad::Exp::PowerOp >, double, false, true, Sacado::Fad::Exp::ExprSpecDefault, Sacado::Fad::Exp::PowerImpl::Scalar>’ and ‘const dealii::Tensor<2, 2, Sacado::Fad::Exp::GeneralFad > >’)  294 |   return std::pow(determinant(F), -1.0 / dim) * F;  |  ~^~~ And some warnings: /home/michael/dealii-9.2.0/examples/code-gallery/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc: In instantiation of ‘void Cook_Membrane::Solid::system_setup() [with int dim = 2; NumberType = double]’:/home/michael/dealii-9.2.0/examples/code-gallery/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc:1005:5:   required from ‘void Cook_Membrane::Solid::run() [with int dim = 2; NumberType = double]’/home/michael/dealii-9.2.0/examples/code-gallery/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc:2263:24:   required from here/home/michael/dealii-9.2.0/examples/code-gallery/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc:1136:35: warning: ‘void dealii::DoFTools::count_dofs_per_block(const DoFHandlerType&, std::vector&, const std::vector&) [with DoFHandlerType = dealii::DoFHandler<2, 2>]’ is deprecated [-Wdeprecated-declarations]1136 | DoFTools::count_dofs_per_block(dof_handler_ref, dofs_per_block,  | ~~^1137 |    block_compone

RE: [deal.II] Matrix Multiplication

[Michael Li]

Hi Jean-Paul, I revisited the code gallery example “The deal.II Library: The 'Quasi-Static Finite-Strain Compressible Elasticity' code gallery program (dealii.org)”. It seems I just need to modify two functions, i.e., getJc(det_F, b_bar) and get_tau(det_F, b_bar) if I want to implement the St. Vienant Kirchhoff material. Is that correct? I compiled cook_membrane.cc a couple of days ago, but for no reason, it fails giving the following errors: /home/michael/dealii-9.2.0/examples/code-gallery/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc:1815:5:   required from here/home/michael/share/trilinos/include/Sacado_Fad_Exp_Ops.hpp:775:1: error: no type named ‘type’ in ‘struct Sacado::mpl::disable_if, Sacado::Fad::Exp::MultiplicationOp >, double, false, true, Sacado::Fad::Exp::ExprSpecDefault, Sacado::Fad::Exp::PowerImpl::Scalar>, double, false, true, Sacado::Fad::Exp::ExprSpecDefault> >’ /usr/local/include/deal.II/physics/elasticity/kinematics.h:294:47: error: no match for ‘operator*’ (operand types are ‘Sacado::Fad::Exp::PowerOp >, double, false, true, Sacado::Fad::Exp::ExprSpecDefault, Sacado::Fad::Exp::PowerImpl::Scalar>’ and ‘const dealii::Tensor<2, 2, Sacado::Fad::Exp::GeneralFad > >’)  294 |   return std::pow(determinant(F), -1.0 / dim) * F;  |  ~^~~ And some warnings: /home/michael/dealii-9.2.0/examples/code-gallery/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc: In instantiation of ‘void Cook_Membrane::Solid::system_setup() [with int dim = 2; NumberType = double]’:/home/michael/dealii-9.2.0/examples/code-gallery/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc:1005:5:   required from ‘void Cook_Membrane::Solid::run() [with int dim = 2; NumberType = double]’/home/michael/dealii-9.2.0/examples/code-gallery/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc:2263:24:   required from here/home/michael/dealii-9.2.0/examples/code-gallery/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc:1136:35: warning: ‘void dealii::DoFTools::count_dofs_per_block(const DoFHandlerType&, std::vector&, const std::vector&) [with DoFHandlerType = dealii::DoFHandler<2, 2>]’ is deprecated [-Wdeprecated-declarations] 1136 | DoFTools::count_dofs_per_block(dof_handler_ref, dofs_per_block,  | ~~^ 1137 |    block_component);  |      Thanks,Michael From: Jean-Paul PelteretSent: Thursday, July 22, 2021 11:03 AMTo: dealii@googlegroups.comSubject: Re: [deal.II] Matrix Multiplication Dear Michael, Although this BDB form is, as you say, ubiquitous in the literature, it is a consequence to using global scalar-valued shape functions everywhere as opposed to vector-valued shape functions that we support. So it falls out of the way that the finite element Ansatz is defined for a vector field. It’s definitely possible to write a deal.II program in this BDB form (a former colleague did this), but its completely not straightforward. I’m afraid to say that you’d have to work out all of the translation from the tensor-based, per DoF shape function operators to this global cell operator yourself.  Have you looked at the finite strain solid mechanics tutorials and code gallery examples? They are really implementing the same thing, only with the native or “natural" deal.II constructs due to the use of vector-valued shape functions. It's easy enough to translate those tutorials to some other parameterisation, if need be. Best,Jean-PaulOn 22. Jul 2021, at 17:47, Michael Li <lianxi...@gmail.com> wrote: Dear Jean-Paul,Thank you for your information. It seems tensor is an elegant way to do the computation. However I don’t know how to translate the following formulation with tensor notation. The Voigt notation is used ubiquitous in the FEM formulation of solid mechanics. It is straightforward to use matrix vector operation as demonstrated in many references. In the following equation, the B matrix (displacement-strain matrix with geometric nonlinearity) is not square, so it looks like Bcannot be represented by a tensor. If I use Voigt notation, then it seems I can only choose FullMatrix class which looks clumsy and not fully exploiting deal.II. I suppose it can be implemented in deal.II using tensor notation but just have no idea how to do that.  <991EE976D38946F2A751BB07719E062B.png> Thanks,Michael  From: Jean-Paul PelteretSent: Wednesday, July 21, 2021 10:45 PMTo: dealii@googlegroups.comSubject: Re: [deal.II] Matrix Multiplication Dear Michael, The classes that better serve your purpose are the Tensor and SymmetricTensor classes:https://dealii.org/current/doxygen/deal.II/classTensor.htmlhttps://dealii.org/current/doxygen/deal.II/classSymmetricTensor.htmlThese classes have a 

RE: [deal.II] Questions on Step-18

[Michael Li]

t; wrote: Hi Jean-Paul and All, I would like to discuss the formulas in step-18 further.   In the code, the angles are computed as follows. For 2D, const double curl = (grad_u[1][0] - grad_u[0][1]);    const double angle = std::atan(curl); For 3D,const double tan_angle = std::sqrt(curl * curl);const double angle = std::atan(tan_angle); I don’t really understand why we need to do the atan operation to get the angle. It seems to me that half of the curl of displacement itself is the rotation angle.  Moreover, for 2D, do we need to do the same as 3D   tan_angle = std::sqrt(curl * curl)? I mean we first calculate the magnitude of the curl then calculate the angle. I’m concerned about the sign of the curl in 2D angle = std::atan(curl) (we need the magnitude of the curl). I wonder what references I can get from these formulas. I hope I can be fully convinced of the stuff.  To be specific, I want to make sure where I should put the factor of ½ to make the patch. I.e., which is correct in the following?  1.  const double tan_angle = std::sqrt(0.5 * curl * 0.5 * curl);or 2.  const double tan_angle = 0.5 * std::sqrt(curl * curl);  In my opinion, the code should be as follows.   For 2D,const double curl = (grad_u[1][0] - grad_u[0][1]); const double angle = 1.0 / 2.0 * std::sqrt(curl * curl); // Half of the magnitude of curl is the rotation angle. For 3D,const double angle = 1.0 / 2.0 * std::sqrt(curl * curl);const double angle = std::atan(tan_angle);  Any comments will be of great help. Best,Michael On Mon, Jun 28, 2021 at 2:21 AM Andrew McBride <mcbride.and...@gmail.com> wrote:My view here is that the problem is quasi-static and hence time serves to order events. Hence a displacement increment can be viewed as equivalent to the velocity.  For example, let’s fix the number of time-steps to 10. If the total time is 1s or 10s we should get the same results (assuming the timing of the application of loading is scaled based on the assumption that time simply order events). Best,AndrewOn 27 Jun 2021, at 02:19, Michael Li <lianxi...@gmail.com> wrote: Hi Jean-Paul,Thank you for the nice explanation and information! It cleared out my concern and helped me understand the related concepts of vorticity and rotation. Vorticity (1/2 curl of velocity) stands for the rotation rate (angular velocity) but the infinitesimal rotation tensor gives the rotation angle (curl of displacement). Do correct me if I’m wrong. Now I’m excited to learn how to make my first patch. - Michael From: Jean-Paul PelteretSent: Saturday, June 26, 2021 2:57 PMTo: dealii@googlegroups.comSubject: Re: [deal.II] Questions on Step-18 Hi Michael, What’s written in the implementation suggests that computing the curl of the displacement (increments) is the right thing to do in this instance: Nevertheless, if the material was prestressed in a certain direction, then this direction will be rotated along with the material. To this end, we have to define a rotation matrix R(Δun) that describes, in each point the rotation due to the displacement increments. It is not hard to see that the actual dependence of R on Δun can only be through the curl of the displacement, rather than the displacement itself or its full gradient (as mentioned above, the constant components of the increment describe translations, its divergence the dilational modes, and the curl the rotational modes).  I think that this is because the displacement field for solids is analogous to the velocity field for fluids, for which that equation in the Wiki link was expressly written (Compare incompressible linear elasticity to Stokes flow — the equations have the same form). Furthermore, having dug through my continuum mechanics notes, I now see that this specific rotation matrix is called the “infinitesimal rotation tensor” and is defined as the antisymmetric part of \grad u. Since step-18 is dealing with incremental updates, the incremental form of the rotation tensor is computed. Sorry for the confusion. So, to my understanding it’s just the factor of 1/2 that needs to be added. Best,Jean-Paul On 26. Jun 2021, at 16:14, Michael Lee <lianxi...@gmail.com> wrote: I would be happy to do the comparison and make the fix. But before doing that, I want to make sure I understand the formula correctly. When we calculate the rotation matrix, it seems that the du (incremental displacement) is used.  // Then initialize the FEValues object on the present  // cell, and extract the gradients of the displacement at the  // quadrature points for later computation of the strains  fe_values.reinit(cell);  fe_values.get_function_gradients(incremental_displacement,   displacement_increment_grads);Should we use the velocity (), since ? Can you check this as well?Best,Michael   On Sat, Jun 26, 2021 at 1:48 AM Jean-Paul Pelteret <jppelte...@gmail.com> wrote:Following Andrew’s expla

RE: [deal.II] Questions on Step-18

[Michael Li]

Hi Jean-Paul,Thank you for the nice explanation and information! It cleared out my concern and helped me understand the related concepts of vorticity and rotation. Vorticity (1/2 curl of velocity) stands for the rotation rate (angular velocity) but the infinitesimal rotation tensor gives the rotation angle (curl of displacement). Do correct me if I’m wrong. Now I’m excited to learn how to make my first patch. - Michael From: Jean-Paul PelteretSent: Saturday, June 26, 2021 2:57 PMTo: dealii@googlegroups.comSubject: Re: [deal.II] Questions on Step-18 Hi Michael, What’s written in the implementation suggests that computing the curl of the displacement (increments) is the right thing to do in this instance: Nevertheless, if the material was prestressed in a certain direction, then this direction will be rotated along with the material. To this end, we have to define a rotation matrix R(Δun) that describes, in each point the rotation due to the displacement increments. It is not hard to see that the actual dependence of R on Δun can only be through the curl of the displacement, rather than the displacement itself or its full gradient (as mentioned above, the constant components of the increment describe translations, its divergence the dilational modes, and the curl the rotational modes).  I think that this is because the displacement field for solids is analogous to the velocity field for fluids, for which that equation in the Wiki link was expressly written (Compare incompressible linear elasticity to Stokes flow — the equations have the same form). Furthermore, having dug through my continuum mechanics notes, I now see that this specific rotation matrix is called the “infinitesimal rotation tensor” and is defined as the antisymmetric part of \grad u. Since step-18 is dealing with incremental updates, the incremental form of the rotation tensor is computed. Sorry for the confusion. So, to my understanding it’s just the factor of 1/2 that needs to be added. Best,Jean-PaulOn 26. Jun 2021, at 16:14, Michael Lee <lianxi...@gmail.com> wrote: I would be happy to do the comparison and make the fix. But before doing that, I want to make sure I understand the formula correctly. When we calculate the rotation matrix, it seems that the du (incremental displacement) is used.  // Then initialize the FEValues object on the present  // cell, and extract the gradients of the displacement at the  // quadrature points for later computation of the strains  fe_values.reinit(cell);  fe_values.get_function_gradients(incremental_displacement,   displacement_increment_grads);Should we use the velocity (), since ? Can you check this as well?Best,Michael   On Sat, Jun 26, 2021 at 1:48 AM Jean-Paul Pelteret <jppelte...@gmail.com> wrote:Following Andrew’s explanation, I suppose that this is relation for which we’re lacking the factor of 1/2, right? https://en.wikipedia.org/wiki/Angular_velocity#Angular_velocity_as_a_vector_field If so, then maybe we should link to this equation in the tutorial documentation too. If this is wrong in the deal.II code, would you be interested in writing a patch to correct this? It should be an easy fix, and a good first patch to contribute to the library! I concur — it would be great if you’d be willing to write a patch that fixes this! It would be interesting to see the effect on the results of the tutorial. Best,Jean-Paul On 26. Jun 2021, at 05:58, Wolfgang Bangerth <bange...@colostate.edu> wrote: On 6/24/21 6:32 PM, Michael Li wrote:Andrew, thanks for confirming that. The missing 1/2 does not affect the demonstration of functionalities of deal.II but it may change the results.If this is wrong in the deal.II code, would you be interested in writing a patch to correct this? It should be an easy fix, and a good first patch to contribute to the library!If you're curious about ho to write a patch, take a look at https://github.com/dealii/dealii/wiki/ContributingBestW.  On Sat, Jun 26, 2021 at 1:48 AM Jean-Paul Pelteret <jppelte...@gmail.com> wrote:Following Andrew’s explanation, I suppose that this is relation for which we’re lacking the factor of 1/2, right? https://en.wikipedia.org/wiki/Angular_velocity#Angular_velocity_as_a_vector_field If so, then maybe we should link to this equation in the tutorial documentation too. If this is wrong in the deal.II code, would you be interested in writing a patch to correct this? It should be an easy fix, and a good first patch to contribute to the library! I concur — it would be great if you’d be willing to write a patch that fixes this! It would be interesting to see the effect on the results of the tutorial. Best,Jean-Paul On 26. Jun 2021, at 05:58, Wolfgang Bangerth <bange...@colostate.edu> wrote: On 6/24/21 6:32 PM, Michael Li wrote:Andrew, thanks for confirming that. The missing 1/2 does not affect the demonstration of functionalities of deal.II bu

RE: [deal.II] Questions on Step-18

[Michael Li]

Hi Jean-Paul,Thank you so much for the deliberate suggestions! It seems the code-gallery example is a good starting point for me to implement St-Venant Kirchhoff hyperelastic model with the help of “push-forward” transformations. I need to dig into this example to understand how the geometric stiffness tangent is embedded for free. It’s good to know a tutorial that uses one-filed elasticity will be available. I believe that will help a lot for beginners. AD is an interesting tool for nonlinear FEM. In  an early discussion, Alex shared his similar codes using AD at Integrated material and spatial traction forces on boundary not equal (google.com). But since I’m a deal.II beginner and not familiar with AD, I have not studied his code further. So I am going to study your code-gallery example first as it is well documented. I may ask for your further help and will let you know if I get the code roll. Thank you again for your wonderful help! Best,Michael  From: Jean-Paul PelteretSent: Thursday, June 24, 2021 11:43 PMTo: dealii@googlegroups.comSubject: Re: [deal.II] Questions on Step-18 HI Michael. To be honest, I have a hard time to understand the three-filed formulation in Step 44.  That’s a fair comment, and a known deficit in our current tutorial list. I hope that by the end of the year we will have a tutorial that uses one-field elasticity, and would be a much better starting point for you and others (WIP: https://github.com/dealii/dealii/pull/10394). However, I’m still interested in implementing the pure geometrical large finite deformation model. I believe that will make things simpler and help me understand more complex nonlinear models. Meanwhile, I try to compare the results with a reference which only consider the geometrical nonlinearity. It sounds to me like you’re looking to implement finite strain elasticity with a St Venant-Kirchhoff constitutive law. That is, as far as I understand, the simplest constitutive law for finite strain elasticity.https://en.wikipedia.org/wiki/Hyperelastic_material#Saint_Venant%E2%80%93Kirchhoff_modelIn this model, the elastic tangent looks like that for linear elasticity, but the strain tensor of interest is no longer the linear/small strain tensor, but rather the (nonlinear) Green-Lagrange strain tensor. (Nevertheless, even though the material law is linear, you still have the geometric tangent contribution that appears after linearisation of the residual, which comes from the linearisation of the now nonlinearly deformation-dependent test function “de” in  \int de : \sigma dv .) If this is indeed what you want, then it is possible to swap out the Neo-Hookean constitutive law that is implemented in the code-gallery example with the St-Venant Kirchhoff constitutive law. You’d get “for free” the geometric stiffness contribution to the tangent, and would just have to worry about the material tangent and definition of the stress. I think that the only complication is the default parameterisation for the constitutive law — it was parameterised in terms of the left Cauchy-Green tensor as the whole problem is set up with the Kirchhoff stress (that’s the Cauchy stress, but per unit reference volume). Nevertheless, with the help of these “push-forward” transformations you would be able to use whatever parameterisation you’d like for the constitutive law (e.g. in terms of the Green-Lagrange strain tensor) and transform the resulting stress and material tangent to the correct (i.e. the spatial) configuration as is required for that implementation of the weak forms. The important thing to remember is that, despite all of the transformations involved, you’re still solving the same BVP with the same constitutive law irrespective of which choice of stress-strain (energetic) conjugate pairs you use, and irrespective of whether its expressed in the total or updated Lagrangian form. But anyway, this is just a suggestion. Irrespective of whether or not you choose to go that route, I wish you luck in implementing what you’re needing to. Best,Jean-PaulOn 25. Jun 2021, at 02:32, Michael Li <lianxi...@gmail.com> wrote: Andrew, thanks for confirming that. The missing 1/2 does not affect the demonstration of functionalities of deal.II but it may change the results. Jean-Paul, thanks for commenting on my second question. I want to study the pure geometrical nonlinear elasticity (large deformation with linear material). It should be the simplest nonlinear model in elasticity. Step 18 looks like a good start as an updated Lagrangian formulation; but it does not include the nonlinear part of the finite strain. I spent some on Step 44 because “ Quasi-Static Finite-Strain Compressible Elasticity” also mentions it is based on Step 44. To be honest, I have a hard time to understand the three-filed formulation in Step 44. So maybe I should get into “Quasi-Static … Elasticity” first which uses classical one-field formulation. However, I’m still interested in implementing the pure geometrical

RE: [deal.II] Questions on Step-18

[Michael Li]

Andrew, thanks for confirming that. The missing 1/2 does not affect the demonstration of functionalities of deal.II but it may change the results. Jean-Paul, thanks for commenting on my second question. I want to study the pure geometrical nonlinear elasticity (large deformation with linear material). It should be the simplest nonlinear model in elasticity. Step 18 looks like a good start as an updated Lagrangian formulation; but it does not include the nonlinear part of the finite strain. I spent some on Step 44 because “ Quasi-Static Finite-Strain Compressible Elasticity” also mentions it is based on Step 44. To be honest, I have a hard time to understand the three-filed formulation in Step 44. So maybe I should get into “Quasi-Static … Elasticity” first which uses classical one-field formulation. However, I’m still interested in implementing the pure geometrical large finite deformation model. I believe that will make things simpler and help me understand more complex nonlinear models. Meanwhile, I try to compare the results with a reference which only consider the geometrical nonlinearity. -Michael   From: Andrew McBrideSent: Thursday, June 24, 2021 2:23 PMTo: deal.II User GroupSubject: Re: [deal.II] Questions on Step-18 Hi both, Jean-Paul has addressed the second point nicely. On the first point, I think there is a 1/2 missing. The curl of the velocity gradient is the vorticity which is twice the angular velocity - hence I think you need a 1/2. Happy to be corrected on this. Best,AndrewOn 24 Jun 2021, at 21:03, Jean-Paul Pelteret  wrote: Hi Michael, I cannot comment on the first question, but might be able to assist a bit with the second. But may I first ask, what precisely are you trying to achieve with this extension?  As interesting as it is, in the past I had found step-18 to deviate too significantly from the “classical” approach to elasticity to be a natural candidate extend towards finite strain elasticity, for example (this is kind-of implied by the caveat that you partially quoted). I’ve been looking at some of my textbooks (e.g. reviewing the topic of the updated Lagrangian formulation in Holzapfel’s “Nonlinear solid mechanics” and Wrigger’s “Nonlinear finite element methods”) to try to answer (2), but cannot trivially reconcile the two approaches. I think that there’s a bit too much going on to be able to correctly deduce by eye what the requisite changes are. I think that you’d need to reformulate the balance laws and consider their implications for the implemented weak forms — in particular, I think that you’d be missing a contribution that (for finite deformation) looks like the geometric stiffness, but there could be further differences that extend from the overall approach taken to the problem.  If you’re interested in an examples that are more closely aligned with what you might see in the literature, then you can take a look at step-44 or the “Quasi-Static Finite-Strain Compressible Elasticity” code-gallery example, which is effectively step-44 reduced to the one-field total Lagrangian formulation that you’d find in many standard textbooks. It would be easy enough to rework this to use the updated Lagrangian approach, if that it what you desire. I hope that this helps you a little. Best,Jean-Paul On 22. Jun 2021, at 15:24, Michael Lee  wrote: Hello,I have two questions when studying Step-18.  1) Should there be a factor 1/2 when calculating the rotation matrix angle (code in tutorial: angle = std::atan(curl) ) ? This is a mathematical problem. Can anyone give some material for the formula to clear out my doubt? 2) In the introduction, it says "we will consider a model in which we produce a small deformation, deform the physical coordinates of the body by this deformation, and then consider the next loading step again as a linear problem. This isn't consistent, since the assumption of linearity implies that deformations are infinitesimal and so moving around the vertices of our mesh by a finite amount before solving the next linear problem is an inconsistent approach." My question is how to make it consistent. Can I just add the nonlinear part of the strain tensor like the following (of course in both get_strain functions)? I just want to consider the geometrical nonlinearity.   template   inline SymmetricTensor<2, dim>  get_strain(const std::vector> )  {Assert(grad.size() == dim, ExcInternalError()); SymmetricTensor<2, dim> strain;for (unsigned int i = 0; i < dim; ++i)  strain[i][i] = grad[i][i]; for (unsigned int i = 0; i < dim; ++i)  for (unsigned int j = i + 1; j < dim; ++j)strain[i][j] = (grad[i][j] + grad[j][i] +       // linear part                              grad[i][0] * grad[j][0] +    //nonlinear part                              grad[i][1] * grad[j][1] +                                 grad[i][2] * grad[j][2]) / 2; return strain;  } Thanks,Michael -- The deal.II project is located at 

RE: [deal.II] Different get_function_gradients

[Michael Li]

Hi Jean-Paul,Thank you so much for your explanation. That totally makes sense. These different functions exist for implementing  different weak forms easily. Though I can get everything using (1), but I would prefer using the others when I try to extract components of the gradient.  Best,Michael From: Jean-Paul PelteretSent: Thursday, June 24, 2021 2:13 PMTo: dealii@googlegroups.comSubject: Re: [deal.II] Different get_function_gradients Hi Michael,Does the (1) do all the work that the last two do? I.e., (1) gives the whole chunk of gradient, and the last two extract the scalar and vector parts from the chunk, respectively.They are very closely related. (1) has no knowledge on the nature of the field(s) for which the gradient is being computed, so it just returns the gradient of all individual components of the solution field at once, and it's up to you to interpret that however you need. For (2) and (3), one is interpreting the solution field (or a subset of its components) to have a certain nature, and so more context can be given to the gradients that are returned. In particular, the gradient of a scalar field is a rank-1 tensor, and the gradient of a vector field would be a rank-2 tensor. So the return types are in fact different, but essentially encode the some/all information that is given by (1).  Does that make sense? The great thing about using the FEValuesExtractors/FEValuesViews concept is that you are able to achieve a near 1-1 match of the written and programmed weak form, which makes implementing things a lot easier. So in most of the more recent tutorials we use them, and I’d say that more often than not we’d advocate their use. Best,Jean-PaulOn 24. Jun 2021, at 18:20, Michael Lee  wrote: After reading the documentation on Handling vector valued problems (The deal.II Library: Handling vector valued problems (dealii.org)), I have the following question.What are the differences between the three functions?(1) FEValuesBase::get_function_gradients (2) FEValuesViews::Scalar::get_function_gradients (3) FEValuesViews::Vector::get_function_gradients  Does the (1) do all the work that the last two do? I.e., (1) gives the whole chunk of gradient, and the last two extract the scalar and vector parts from the chunk, respectively.  -- 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/89ec12a9-9553-4ee2-9bf0-8900517587c3n%40googlegroups.com. -- 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/97AA7B43-5CF1-4ECE-9238-272250916C76%40gmail.com. 



-- 
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/A647AB72-A22E-46CB-9AF9-06FCEA7C36C4%40hxcore.ol.

RE: [deal.II] Re: Integrated material and spatial traction forces on boundary not equal

[Michael Li]

Thanks, W. for confirming that! I reinstalled deal.II after installing Trillinos. Step 31 runs well.   From: Wolfgang BangerthSent: Wednesday, June 16, 2021 4:10 PMTo: dealii@googlegroups.comSubject: Re: [deal.II] Re: Integrated material and spatial traction forces on boundary not equal On 6/16/21 4:07 PM, Michael Lee wrote:> > Is there a way to let deal.II know Trillinos is installed? Do I have to > reinstall deal.ii?> Or, can I compile the code without using Trillinos with some modification? No, you have to re-install deal.II so that at compile time, it knows that Trilinos exists. Best  W. -- Wolfgang Bangerth  email: bange...@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 a topic in the Google Groups "deal.II User Group" group.To unsubscribe from this topic, visit https://groups.google.com/d/topic/dealii/K-lMxbtZUdQ/unsubscribe.To unsubscribe from this group and all its topics, send an email to dealii+unsubscr...@googlegroups.com.To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/06a30833-0914-ed13-b4e1-dfd9c2c3b72d%40colostate.edu. 



-- 
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/F7A11E7B-B4F1-4CEA-B54F-01D974A17F04%40hxcore.ol.

RE: [deal.II] How to define boundary conditions on edges and vertices in 3D?

[Michael Li]

Thanks for confirming that. It saves me from keeping trying to define boundary conditions by setting boundary id on edges or vertices. I was too excited to think it is possible to do that when reading the code line(j)->set_boundary_id(2) in grid_generator.cc. Maybe it is just to distinct the edges with different boundary_ids for other use. There are still many scenarios that we want to fix some points in the domain. I wonder if we can make use of the AffineConstraints class to deal with Dirichlet boundary conditions though it is not possible for Neumman BCs. Best,Michael From: Wolfgang BangerthSent: Thursday, June 10, 2021 10:41 AMTo: dealii@googlegroups.comSubject: Re: [deal.II] How to define boundary conditions on edges and vertices in 3D? On 6/10/21 8:31 AM, Michael Lee wrote:> 1) thin plate under pressure with *four edges simply supported*.> 2) thin plate under pressure with *four corner points simply supported*.> 2021-06-10_082227.png> For the first one, I tried the following code, but it seems boundary condition > only works with face not on edge. That is correct. Mathematically, you can only prescribe boundary conditions on portions of the boundary (i.e., on *faces*) but not on edges in 3d or vertices in 2d/3d. What is easy to do in deal.II matches what is mathematically possible. Now, you will ask why other software packages allow you to do what is not easy in deal.II. The answer is that other packages generally use a fixed mesh, and on a fixed mesh, prescribing boundary conditions on a vertex is roughly equivalent to prescribing boundary conditions on the adjacent faces. That's ok if you keep the mesh fixed, but it means that you are changing the kind of boundary conditions every time you change the mesh if you use software that uses adaptive mesh refinement. So my suggestion would be to think about the physical situation you want to model, and translate that into prescribing boundary conditions on faces instead of edges/vertices. Best  W. -- Wolfgang Bangerth  email: bange...@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/9e09474b-1e3f-7fbf-7527-01f6cdd256ef%40colostate.edu. 



-- 
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/89FC6FE0-6F24-4447-B995-03F9B36D0E84%40hxcore.ol.

RE: [deal.II] Re: Integrated material and spatial traction forces on boundary not equal

[Michael Li]

Alex,  Thank you for sharing your codes. I have some compiling errors relating to “EnergyFunctional’: solve_ring_nonlinear.cpp:520:43: error: ‘EnergyFunctional’ in namespace ‘dealii::Differentiation::AD’ does not name a template type  520 | using ADHelper = Differentiation::AD::EnergyFunctional<  |   ^~~~Can you help me with that? For the large error, I noticed you’re using linear element. I encountered the same large error when comparing it with Abaqus with FE_Q(1). But the error came down with less grids when I used higher finite element FE_Q(2). I remember the deflection of beam is a cubic function of coordinate. You may try that see if it improves. Best,Michael From: Alex CumberworthSent: Wednesday, June 9, 2021 9:54 AMTo: deal.II User GroupSubject: [deal.II] Re: Integrated material and spatial traction forces on boundary not equal Hello, I have attached the most recent version of my code here. I have tried to make setting boundary conditions in the parameter file more convenient for myself; you can set boundary domains, boundary conditions that use these domains, and stages if the displacement is large. However, there are not many comments, so you may want to just remove this part for your own purposes. However, I have been a bit surprised in my comparisons at how fine a mesh is required to achieve convergence with the beam theory result. I am now using a beam that is 1 x 2 x 20, and using the subdivided rectangle helper function, I set the number of subdivisions to be 5 and 10 for the width and height, respectively. I then varied the number of subdivision in the length between 10 and 1000. The beam theory result is that the shear force has a magnitude of 0.001 for a displacement on the right side of 0.1. Even at 1000 subdivisions, the FEM result is 0.00113 (from 0.00129 at 500). The system has 200 000 degrees of freedom, and the result is still off by 13%. Is it expected that even to calculate a shear force in this simple problem that such a large number of degrees of freedom are required?  Best,Alex On Monday, June 7, 2021 at 3:39:44 p.m. UTC+2 lian...@gmail.com wrote:Hi Alex,I'm learning deal.ii and trying do the similar verification. If it is possible for you to share the code with me? Thank you!MichaelOn Tuesday, May 11, 2021 at 4:46:55 AM UTC-6 alexanderc...@gmail.com wrote:Hello, As a test to validate my code, I am solving the equations for geometrically nonlinear elasticity (the Saint Venant-Kirchhoff model) for a beam with a small displacement boundary condition on the right end such that I can compare with Euler-Bernoulli beam theory. I can compare both the displacement and the shear force between the FEM solution and the beam theory solution. In my FEM integration, I output the normal and shear forces for both sides of the beam in both the material and spatial reference. The left and right sides are balanced, as expected, but the spatial and material forces are not quite equal. Shouldn't it be the case that spatial and material force is the same? Here are the outputted forces for the right side Right boundary material normal force: 0.0694169Right boundary spatial normal force: 0.0724468Right boundary material shear force: 0.152057Right boundary spatial shear force: 0.152864 Further, beam theory gives a shear force with a magnitude of exactly 0.2. If I make the displacement smaller the FEM and beam theory shear forces do not converge. Is it expected for them to converge? Below is the deformed system with the stress vectors on the faces included. The black grid is the deformed FEM solution, while the solid red is the beam theory solution.If there is an issue, I would guess it would be in the integration. In converting the material normal vector to the spatial reference, I first only applied the inverse transpose of the deformation gradient, and did not multiply by the determinant until calculating the force vector. I did this so that I can get the unit normal spatial vectors to add up and later average so that I have an average normal vector for the whole boundary face to calculate the normal and shear force vectors. I have pasted the function in below: template void SolveRing::integrate_over_boundaries() {    QGauss quadrature_formula(fe.degree + 1);    FEFaceValues fe_face_values(    fe,    quadrature_formula,    update_values | update_gradients | update_quadrature_points |    update_JxW_values | update_normal_vectors);    std::vector> material_force(2);    std::vector> spatial_force(2);    std::vector> ave_material_normal(2);    std::vector> ave_spatial_normal(2);    const FEValuesExtractors::Vector displacements(0);    for (const auto& cell: dof_handler.active_cell_iterators()) {    for (const auto face_i: GeometryInfo::face_indices()) {    const unsigned int boundary_id {cell->face(face_i)->boundary_id()};    if (not(boundary_id == 1 or boundary_id == 2)) {