# [petsc-users] Elasticity tensor in ex52

Miguel Angel Salazar de Troya salazardetroya at gmail.com
Mon Apr 21 10:27:55 CDT 2014

Thank you.

On Mon, Apr 21, 2014 at 10:24 AM, Matthew Knepley <knepley at gmail.com> wrote:

> On Mon, Apr 21, 2014 at 10:23 AM, Miguel Angel Salazar de Troya <
> salazardetroya at gmail.com> wrote:
>
>> Thank you. Now it makes sense. Just to confirm though, back to the weak
>> form I wrote above in index notation:
>>
>> phi_{i,j} C_{ijkl} u_{k,l}
>>
>> The g3 indices are equivalent to the indices of this weak form in this
>> manner:
>>
>> ic = i
>> jc = k
>> id = j
>> jd = l
>>
>> Is this correct? This is what I understand when you talk about the
>> components and the derivatives of the components.
>>
>
> Yep, that is it.
>
>   Thanks,
>
>      Matt
>
>
>> On Mon, Apr 21, 2014 at 9:48 AM, Matthew Knepley <knepley at gmail.com>wrote:
>>
>>> On Sun, Apr 20, 2014 at 4:17 PM, Miguel Angel Salazar de Troya <
>>> salazardetroya at gmail.com> wrote:
>>>
>>>> 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.
>>>>
>>>
>>> I made an error explaining g3, which is indexed
>>>
>>>   g3[ic, jc, id, jd]
>>>
>>> I thought this might be better since, it is composed of dim x dim
>>> blocks. I am not opposed to changing
>>> this if there is evidence that another thing is better.
>>>
>>>    Matt
>>>
>>>
>>>>
>>>> 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?
>>>>>>>>
>>>>>>>> Miguel
>>>>>>>>
>>>>>>>> --
>>>>>>>> *Miguel Angel Salazar de Troya*
>>>>>>>> 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
>>>>>>> -- Norbert Wiener
>>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>> --
>>>>>> *Miguel Angel Salazar de Troya*
>>>>>> 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
>>>>> -- Norbert Wiener
>>>>>
>>>>
>>>
>>>
>>> --
>>> What most experimenters take for granted before they begin their
>>> experiments is infinitely more interesting than any results to which their
>>> -- Norbert Wiener
>>>
>>
>>
>>
>> --
>> *Miguel Angel Salazar de Troya*
>> 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
> -- Norbert Wiener
>

--
*Miguel Angel Salazar de Troya*