Question about how to get the Stiffness matrix partial derivative with respect to density

[Lance Zhang]

Hello ,

Thanks for your time to see this question.

The given conditions are stiffness matrix,destination function and x from Ax=b,x is movement ,b is force,A is global stiffness matrix.

J(x)=||x||norm-1

There is one equation:

dJ(x,θ)/dθ=∂J(x,θ)/∂θ+∂J(x,θ)/∂x *∂x/∂θ=∂J(x,θ)/∂θ+∂J(x,θ)/∂x*[A^-1(θ))(-∂A(θ)/∂θ)x],because b has no relation with density.

I don't know how to get [A^-1(θ))(-∂A(θ)/∂θ)x],A is matrix with nxn,if θ is mx1 vector,∂A(θ)/∂θ is nxnxm matrix. 

But I don't know how to calculate (A^-1(θ))(-∂A(θ)/∂θ),I can simply get the value of A^-1.

I use cook_membrane.cc this file for dJ(x,θ)/dθ.

I found the density was just included from the mu invoked from parameters.prm.

>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

public:
    Material_Compressible_Neo_Hook_One_Field(const double mu,
                                             const double nu)
      :
      kappa((2.0 * mu * (1.0 + nu)) / (3.0 * (1.0 - 2.0 * nu))),
      c_1(mu / 2.0)
    {
      Assert(kappa > 0, ExcInternalError());
    }

    const double kappa;
    const double c_1;

    // Value of the volumetric free energy
    NumberType
    get_Psi_vol(const NumberType &det_F) const
    {
      return (kappa / 4.0) * (det_F*det_F - 1.0 - 2.0*std::log(det_F));
    }

    NumberType
    get_dPsi_vol_dJ(const NumberType &det_F) const
    {
      return (kappa / 2.0) * (det_F - 1.0 / det_F);
    }

    get_d2Psi_vol_dJ2(const NumberType &det_F) const
    {
      return ( (kappa / 2.0) * (1.0 + 1.0 / (det_F * det_F)));
    }    NumberType
    get_Psi_iso(const SymmetricTensor<2,dim,NumberType> &b_bar) const
    {
      return c_1 * (trace(b_bar) - dim);
    }
    SymmetricTensor<2,dim,NumberType>
    get_tau_bar(const SymmetricTensor<2,dim,NumberType> &b_bar) const
    {
      return 2.0 * c_1 * b_bar;
    }

>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

 assemble_system_one_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
                             ScratchData_ASM &scratch,
                             PerTaskData_ASM &data)
    {
      // Due to the C++ specialization rules, we need one more
      // level of indirection in order to define the assembly
      // routine for all different number. The next function call
      // is specialized for each NumberType, but to prevent having
      // to specialize the whole class along with it we have inlined
      // the definition of the other functions that are common to
      // all implementations.
      assemble_system_tangent_residual_one_cell(cell, scratch, data);
      assemble_neumann_contribution_one_cell(cell, scratch, data);
    }

<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<



Unfortunately ,I did not find density vector in this file  to calculate -∂A(θ)/∂θ.

Could anyone provide any hint or suggestions?

Thanks in advance!

Best regards

Lance

[Lance Zhang]

Hello,

I found the information from one website, it could be better  for understanding the equation .

This equation is an adjoint equation.

b is force which has no relation with density.

Thanks in advance!

Best regards

Lance

Lance Zhang <dim...@gmail.com> 于2023年7月24日周一 14:49写道：

--
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/5IN4HUV1UhY/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dealii+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/f67b2b60-c697-4235-bc57-d0db7890dfccn%40googlegroups.com.

[Wolfgang Bangerth]

> The given conditions are stiffness matrix,destination function and x from
> Ax=b,x is movement ,b is force,A is global stiffness matrix.
> J(x)=||x||norm-1

This is a questionable choice because dJ/dx is not defined at x=0.

> There is one equation:
>
> dJ(x,θ)/dθ=∂J(x,θ)/∂θ+∂J(x,θ)/∂x
> *∂x/∂θ=∂J(x,θ)/∂θ+∂J(x,θ)/∂x*[A^-1(θ))(-∂A(θ)/∂θ)x],because b has no relation
> with density.
>
> I don't know how to get [A^-1(θ))(-∂A(θ)/∂θ)x],A is matrix with nxn,if θ is
> mx1 vector,∂A(θ)/∂θ is nxnxm matrix.

∂A(θ)/∂θ)x is an nxm matrix. You want to apply A^{-1} to it. The way you do
that is to call the result
Y = A^{-1} [∂A(θ)/∂θ)x]
and then each column of Y, call it y_k (k=1...m) is computed as
y_k = A^{-1} {[∂A(θ)/∂θ)x]_k}
which requires solving a linear system
A y_k = {[∂A(θ)/∂θ)x]_k}

You never need A^{-1} explicitly.

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Lance Zhang]

Hello Wolfgang,

thanks for your reply.

I will follow your point to see if I could find the solution.

One moire question,how could I get  ∂A(θ)/∂θx,because I did not find any information about density vector like  θ=[θ1,θ2,...,θm],may I know if I have to set density vector value in this finite element?

Code Link:https://www.dealii.org/current/doxygen/code-gallery/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc

mu is invoked from prm file at the position of line 252

Best regards

Lance

[Lance Zhang]

Hello Wolfgang,

The value of  A^-1(θ))(-∂A(θ)/∂θ)x is what I would like to obtain.

Thanks in advance!

Best regards

Lance

[Wolfgang Bangerth]

On 7/24/23 11:24, Lance Zhang wrote:
>
> One moire question,how could I get  ∂A(θ)/∂θx,because I did not find any
> information about density vector like  θ=[θ1,θ2,...,θm],may I know if I have
> to set density vector value in this finite element?

Lance -- I have no idea. I don't know the program, but looking at it, I see no
mention of the word "density" anywhere. Apparently, the person who wrote the
program implemented an equation in which density does not appear. If you want
that the density is part of the formulation, you should first write down what
the specific problem is you want to solve. We can't know what it is you want
to do. You should then compare what you want to do with what the existing
program does, and change it so that it matches to what you want it to do.

[Lance Zhang]

Hello Wolfgang,

thanks for your reply.

Currently,I would like to get the sensitivity by using gradient of objective function with respect to density.

The Objective function is a vector of movement which is calculated from Ax=b.

J=|x|norm1 is the objective function.

dJ(x,Θ)/dΘ is the gradient that we will obtain in the end.

We use adjoint method written in previous emails to calculate dJ(x,Θ)/dΘ.

Best regards

Lance

Wolfgang Bangerth <bang...@colostate.edu> 于 2023年7月24日周一 21:49写道：

--
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/5IN4HUV1UhY/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dealii+un...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/8819fb69-2fc9-a80c-d37d-a45411d0a89b%40colostate.edu.

[blais...@gmail.com]

Dear Lance,

It would be easier if you posted the PDE you wish to solve and it's weak form. This would enable us to help you write the Jacobian for it.

Currently, it is difficult to fully understand what you are trying to achieve.

Best

Bruno