Shape function derivatives

[Seyed Ali Mohseni]

Hi,

I found some issue when constructing the FE B-operator for multiple elements. In a single element it works perfectly.

My approach is the following:

template<int dim>

FullMatrix<double> get_b_operator (const FEValues<dim> &fe_values, const unsigned int dofs_per_cell, const unsigned int q)

{

FullMatrix<double> tmp(dim, GeometryInfo<dim>::vertices_per_cell);

std::vector<DerivativeForm<1, dim, dim> > invJ = fe_values.get_inverse_jacobians();

for (unsigned int i = 0; i < dofs_per_cell; i = i + dim)

{

const unsigned int index = i / dim;

// COMPUTE: B-operator (Remark: This version has to be extended for 3D!)

tmp[0][index] = invJ[q][0][0] * fe_values.shape_grad(i, q)[0] + invJ[q][0][1] * fe_values.shape_grad(i, q)[1];

tmp[1][index] = invJ[q][1][0] * fe_values.shape_grad(i, q)[0] + invJ[q][1][1] * fe_values.shape_grad(i, q)[1];

}

return tmp;

}

Hence, I am computing B = inv(J) * transpose(dN). For a single element I receive the following correct results for dN:

C++ CODE

SHAPE FUNCTION DERIVATIVES (GAUSS POINT 1)

-0.788675  -0.788675  

0.788675  -0.211325  

0.211325  0.211325  

-0.211325  0.788675  

SHAPE FUNCTION DERIVATIVES (GAUSS POINT 2)

-0.788675  -0.211325  

0.788675  -0.788675  

0.211325  0.788675  

-0.211325  0.211325  

SHAPE FUNCTION DERIVATIVES (GAUSS POINT 3)

-0.211325  -0.788675  

0.211325  -0.211325  

0.788675  0.211325  

-0.788675  0.788675  

SHAPE FUNCTION DERIVATIVES (GAUSS POINT 4)

-0.211325  -0.211325  

0.211325  -0.788675  

0.788675  0.788675  

-0.788675  0.211325  

DEAL.II

SHAPE FUNCTION DERIVATIVES

  -0.788675 -0.788675   -0.788675 -0.788675   0.788675 -0.211325   0.788675 -0.211325   -0.211325 0.788675   -0.211325 0.788675   0.211325 0.211325   0.211325 0.211325 

  -0.788675 -0.211325   -0.788675 -0.211325   0.788675 -0.788675   0.788675 -0.788675   -0.211325 0.211325   -0.211325 0.211325   0.211325 0.788675   0.211325 0.788675 

  -0.211325 -0.788675   -0.211325 -0.788675   0.211325 -0.211325   0.211325 -0.211325   -0.788675 0.788675   -0.788675 0.788675   0.788675 0.211325   0.788675 0.211325 

  -0.211325 -0.211325   -0.211325 -0.211325   0.211325 -0.788675   0.211325 -0.788675   -0.788675 0.211325   -0.788675 0.211325   0.788675 0.788675   0.788675 0.788675 

Due to the ordering of the element connectivity (numbering of nodes are different in deal.II) we have to swap column 3 and 4 for comparison.

Now for 4 elements there is something strange when I compute my shape function derivatives within deal.II:

DEAL.II

SHAPE FUNCTION DERIVATIVES

  -1.577350 -1.577350   -1.577350 -1.577350   1.577350 -0.422650   1.577350 -0.422650   -0.422650 1.577350   -0.422650 1.577350   0.422650 0.422650   0.422650 0.422650 

  -1.577350 -0.422650   -1.577350 -0.422650   1.577350 -1.577350   1.577350 -1.577350   -0.422650 0.422650   -0.422650 0.422650   0.422650 1.577350   0.422650 1.577350 

  -0.422650 -1.577350   -0.422650 -1.577350   0.422650 -0.422650   0.422650 -0.422650   -1.577350 1.577350   -1.577350 1.577350   1.577350 0.422650   1.577350 0.422650 

  -0.422650 -0.422650   -0.422650 -0.422650   0.422650 -1.577350   0.422650 -1.577350   -1.577350 0.422650   -1.577350 0.422650   1.577350 1.577350   1.577350 1.577350 

 The results for one element change now by a factor of 2 and with increasing element size it always differs by another factor of 2 again.

I was wondering, if my approach to compute my B-operator was wrong. Doesn't fe_values.shape_grad() give me the shape function derivatives directly? Or do I have to use shape_grad_component? 

Additionally, I checked the symmetric_gradient from step-18 and it contains actually the correct values I need, but somehow in a strange order:

SYMMETRIC GRADIENT (GAUSS POINT 1)

0.000000 0.211325 

0.211325 0.422650

SYMMETRIC GRADIENT (GAUSS POINT 2)

0.000000 0.211325 

0.211325 1.577350

SYMMETRIC GRADIENT (GAUSS POINT 3)

0.000000 0.788675 

0.788675 0.422650 

SYMMETRIC GRADIENT (GAUSS POINT 4)

0.000000 0.788675 

0.788675 1.577350

Kind regards,

S. A. Mohseni

[Wolfgang Bangerth]

On 02/01/2017 07:02 AM, 'Seyed Ali Mohseni' via deal.II User Group wrote:
> My approach is the following:
>
> |
> template<int dim>
> FullMatrix<double> get_b_operator (const FEValues<dim> &fe_values, const
> unsigned int dofs_per_cell, const unsigned int q)
> {
> FullMatrix<double> tmp(dim, GeometryInfo<dim>::vertices_per_cell);
>
> std::vector<DerivativeForm<1, dim, dim> > invJ =
> fe_values.get_inverse_jacobians();
>
> for (unsigned int i = 0; i < dofs_per_cell; i = i + dim)
> {
> const unsigned int index = i / dim;
>
> // COMPUTE: B-operator (Remark: This version has to be extended for 3D!)
> tmp[0][index] = invJ[q][0][0] * fe_values.shape_grad(i, q)[0] +
> invJ[q][0][1] * fe_values.shape_grad(i, q)[1];
> tmp[1][index] = invJ[q][1][0] * fe_values.shape_grad(i, q)[0] +
> invJ[q][1][1] * fe_values.shape_grad(i, q)[1];
> }
>
> return tmp;
> }
> |

fe_values.shape_grad() returns the gradient *on the actual cell*, not on
the reference cell. Consequently, you do not have to multiply it by invJ
again. This is the source of your factor of two.

Best
W.

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

[Seyed Ali Mohseni]

@Prof. Bangerth: Thank you very much :) Finally, I understand how the stiffness matrix in step-18 is constructed. Your shape function derivatives are basically already the B-operator. 

[Wolfgang Bangerth]

On 02/01/2017 09:08 AM, 'Seyed Ali Mohseni' via deal.II User Group wrote:
> @Prof. Bangerth: Thank you very much :) Finally, I understand how the
> stiffness matrix in step-18 is constructed. Your shape function
> derivatives are basically already the B-operator.

Correct. They are rows or columns of that operator.