Re: [deal.II] Shape function derivatives

2017-02-01 Thread Wolfgang Bangerth 

On 02/01/2017 07:02 AM, 'Seyed Ali Mohseni' via deal.II User Group wrote:

My approach is the following:

|
template
FullMatrix get_b_operator (const FEValues _values, const
unsigned int dofs_per_cell, const unsigned int q)
{
FullMatrix tmp(dim, GeometryInfo::vertices_per_cell);

std::vector > 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: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Shape function derivatives

2017-02-01 Thread 'Seyed Ali Mohseni' via deal.II User Group 
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
FullMatrix get_b_operator (const FEValues _values, const 
unsigned int dofs_per_cell, const unsigned int q)
{
FullMatrix tmp(dim, GeometryInfo::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.00 0.211325 
0.211325 0.422650

SYMMETRIC GRADIENT (GAUSS POINT 2)

0.00 0.211325 
0.211325 1.577350

SYMMETRIC GRADIENT (GAUSS POINT 3)

0.00 0.788675 
0.788675 0.422650 

SYMMETRIC GRADIENT (GAUSS POINT 4)

0.00 0.788675 
0.788675 1.577350




Kind regards,
S. A. Mohseni


-- 
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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.