Matrix Multiplication

84 views
Skip to first unread message

Michael Lee

unread,
Jul 21, 2021, 6:39:46 PM7/21/21
to dea...@googlegroups.com
Dear dealii users,

I want to construct a matrix E=1/2(F^TF-1). My approach is as follows.

FullMatrix<double> F(3,3) = 0.0; // deformation gradient
FullMatrix<double> E(3,3) = 0.0; // Green strain
FullMatrix<double> I (IdentityMatrix(3));
F.Tmmult(E, F);
E.add(-1, I);
E.add(2, E);

I wonder if there is an more elegant or easier way to do the matrix construction as I need to do other similar operations.

Best,
Michael

Jean-Paul Pelteret

unread,
Jul 22, 2021, 12:45:31 AM7/22/21
to dea...@googlegroups.com
Dear Michael,

The classes that better serve your purpose are the Tensor and SymmetricTensor classes:
These classes have a number of member and non-member functions that make performing the sort of operations that you listed much easier.
There is also the unit_symmetric_tensor() function that returns the rank-2 identity tensor.

Furthermore, we have a function that will compute the Green-Lagrange strain tensor from a given deformation gradient tensor:

I hope that this helps.

Best,
Jean-Paul

--
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+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/CAEyr2zhNC9RpZVjD%2B4ggHp-PCBCRz9YCtBLSyxxqqeCDE-giSA%40mail.gmail.com.

Michael Li

unread,
Jul 22, 2021, 11:47:41 AM7/22/21
to dea...@googlegroups.com

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 B cannot 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.

 

                        

 

Thanks,

Michael

Jean-Paul Pelteret

unread,
Jul 22, 2021, 1:03:15 PM7/22/21
to dea...@googlegroups.com
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-Paul

On 22. Jul 2021, at 17:47, Michael Li <lian...@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>

Michael Li

unread,
Jul 22, 2021, 5:54:16 PM7/22/21
to dea...@googlegroups.com

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::mpl::is_same<double, double>, Sacado::Fad::Exp::MultiplicationOp<Sacado::Fad::Exp::PowerOp<Sacado::Fad::Exp::GeneralFad<Sacado::Fad::Exp::DynamicStorage<double, double> >, 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<Sacado::Fad::Exp::GeneralFad<Sacado::Fad::Exp::DynamicStorage<double, double> >, double, false, true, Sacado::Fad::Exp::ExprSpecDefault, Sacado::Fad::Exp::PowerImpl::Scalar>’ and ‘const dealii::Tensor<2, 2, Sacado::Fad::Exp::GeneralFad<Sacado::Fad::Exp::DynamicStorage<double, double> > >’)

  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<dim, NumberType>::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<dim, NumberType>::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<unsigned int>&, const std::vector<unsigned int>&) [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);

      |                                    ~~~~~~~~~~~~~~~~

Jean-Paul Pelteret

unread,
Jul 26, 2021, 2:46:34 AM7/26/21
to dea...@googlegroups.com
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.
The 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


Michael Li

unread,
Jul 26, 2021, 12:31:28 PM7/26/21
to dea...@googlegroups.com

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

Reply all
Reply to author
Forward
0 new messages