[petsc-users] Elasticity tensor in ex52

[Miguel Angel Salazar de Troya]

I understand. So if we had a linear elastic material, the weak form would
be something like

phi_{i,j} C_{ijkl} u_{k,l}

where the term C_{ijkl} u_{k,l} would correspond to the term f_1(U,
gradient_U) in equation (1) in your paper I mentioned above, and g_3 would
be C_{ijkl}. Therefore the indices {i,j} would be equivalent to the indices
"ic, id" you mentioned before and "jc" and "jd" would be {k,l}?

For g3[ic, id, jc, jd], transforming the four dimensional array to linear
memory would be like this:

g3[((ic*Ncomp+id)*dim+jc)*dim+jd] = 1.0;

where Ncomp and dim are equal to the problem's spatial dimension.

However, in the code, there are only two loops, to exploit the symmetry of
the fourth order identity tensor:

  for (compI = 0; compI < Ncomp; ++compI) {
    for (d = 0; d < dim; ++d) {
      g3[((compI*Ncomp+compI)*dim+d)*dim+d] = 1.0;
    }
  }

Therefore the tensor entries that are set to one are:

g3[0,0,0,0]
g3[1,1,0,0]
g3[0,0,1,1]
g3[1,1,1,1]

This would be equivalent to the fourth order tensor \delta_{ij}
\delta_{kl}, but I think the one we need is \delta_{ik} \delta_{jl},
because it is the derivative of a matrix with respect itself (or the
derivative of a gradient with respect to itself). This is assuming the
indices of g3 correspond to what I said.

Thanks in advance.

Miguel


On Apr 19, 2014 6:19 PM, "Matthew Knepley" <knepley at gmail.com> wrote:

> On Sat, Apr 19, 2014 at 5:25 PM, Miguel Angel Salazar de Troya <
> salazardetroya at gmail.com> wrote:
>
>> Thanks for your response. Now I understand a bit better your
>> implementation, but I am still confused. Is the g3 function equivalent to
>> the f_{1,1} function in the matrix in equation (3) in this paper:
>> http://arxiv.org/pdf/1309.1204v2.pdf ? If so, why would it have terms
>> that depend on the trial fields? I think this is what is confusing me.
>>
>
> Yes, it is. It has no terms that depend on the trial fields. It is just
> indexed by the number of components in that field. It is
> a continuum thing, which has nothing to do with the discretization. That
> is why it acts pointwise.
>
>    Matt
>
>
>> On Sat, Apr 19, 2014 at 11:35 AM, Matthew Knepley <knepley at gmail.com>wrote:
>>
>>> On Fri, Apr 18, 2014 at 1:23 PM, Miguel Angel Salazar de Troya <
>>> salazardetroya at gmail.com> wrote:
>>>
>>>> Hello everybody.
>>>>
>>>> First, I am taking this example from the petsc-dev version, I am not
>>>> sure if I should have posted this in another mail-list, if so, my
>>>> apologies.
>>>>
>>>> In this example, for the elasticity case, function g3 is built as:
>>>>
>>>> void g3_elas(const PetscScalar u[], const PetscScalar gradU[], const
>>>> PetscScalar a[], const PetscScalar gradA[], const PetscReal x[],
>>>> PetscScalar g3[])
>>>> {
>>>>   const PetscInt dim   = spatialDim;
>>>>   const PetscInt Ncomp = spatialDim;
>>>>   PetscInt       compI, d;
>>>>
>>>>   for (compI = 0; compI < Ncomp; ++compI) {
>>>>     for (d = 0; d < dim; ++d) {
>>>>       g3[((compI*Ncomp+compI)*dim+d)*dim+d] = 1.0;
>>>>     }
>>>>   }
>>>> }
>>>>
>>>> Therefore, a fourth-order tensor is represented as a vector. I was
>>>> checking the indices for different iterator values, and they do not seem to
>>>> match the vectorization that I have in mind. For a two dimensional case,
>>>> the indices for which the value is set as 1 are:
>>>>
>>>> compI = 0 , d = 0      ----->     index = 0
>>>> compI = 0 , d = 1      ----->     index = 3
>>>> compI = 1 , d = 0      ----->     index = 12
>>>> compI = 1 , d = 1      ----->     index = 15
>>>>
>>>> The values for the first and last seem correct to me, but they other
>>>> two are confusing me. I see that this elasticity tensor (which is the
>>>> derivative of the gradient by itself in this case) would be a four by four
>>>> identity matrix in its matrix representation, so the indices in between
>>>> would be 5 and 10 instead of 3 and 12, if we put one column on top of each
>>>> other.
>>>>
>>>
>>> I have read this a few times, but I cannot understand that you are
>>> asking. The simplest thing I can
>>> respond is that we are indexing a row-major array, using the indices:
>>>
>>>   g3[ic, id, jc, jd]
>>>
>>> where ic indexes the components of the trial field, id indexes the
>>> derivative components,
>>> jc indexes the basis field components, and jd its derivative components.
>>>
>>>    Matt
>>>
>>>
>>>> I guess my question is then, how did you vectorize the fourth order
>>>> tensor?
>>>>
>>>> Thanks in advance
>>>> Miguel
>>>>
>>>> --
>>>> *Miguel Angel Salazar de Troya*
>>>> Graduate Research Assistant
>>>> Department of Mechanical Science and Engineering
>>>> University of Illinois at Urbana-Champaign
>>>> (217) 550-2360
>>>> salaza11 at illinois.edu
>>>>
>>>>
>>>
>>>
>>> --
>>> What most experimenters take for granted before they begin their
>>> experiments is infinitely more interesting than any results to which their
>>> experiments lead.
>>> -- Norbert Wiener
>>>
>>
>>
>>
>> --
>> *Miguel Angel Salazar de Troya*
>> Graduate Research Assistant
>> Department of Mechanical Science and Engineering
>> University of Illinois at Urbana-Champaign
>> (217) 550-2360
>> salaza11 at illinois.edu
>>
>>
>
>
> --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which their
> experiments lead.
> -- Norbert Wiener
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20140420/22d5dcfb/attachment.html>