[deal.II] Evaluating the gradient of the solution at the centroid of the element

2016-12-21 Thread cmhamel32
I have gone through some of the tutorials previously and thought it was 
about time I actually tried to implement something on my own in the 
library. 

I am attempting to implement a finite strain code using the F-Bar method. 
Essentially this is a single field formulation which requires the 
calculation the of the deformation gradient at the element center. The 
trick here is though that the element centroid isn't a vertex.

My question is, is there a way to "interpolate" the displacement gradient 
at the element center without explicitly making it a quadrature point or 
will I have to implement my own quadrature formula for this method?

Thank you for any assistance that can be provided.

-- 
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] Re: Postprocessing on Vector Valued Problem

2016-12-21 Thread Jaekwang Kim
You were right

Oh I just got the result! 
I didn't know that dof_handler can automatically matching the sequence of 
the cell. 
Thank you for your time

Jaekwang Kim 

-- 
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] Re: Postprocessing on Vector Valued Problem

2016-12-21 Thread Bruno Turcksin
Jaekwang,

On Wednesday, December 21, 2016 at 4:05:54 PM UTC-5, Jaekwang Kim wrote:
>
>
>  but How can I relate my 'cellwise_shear_rate' with dof handler so that I 
> finally able to see the values in 'vtk' format? 
> for now, it just seems that no more than a list of arrays
>
You mean a list of double not arrays right? It should work fine though, the 
first entry corresponds to the first cell, the second to the second cell, 
etc. Just use data_out.add_vector like usual.

Best,

Bruno 

-- 
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.


Re: [deal.II] deal.II FE structure

2016-12-21 Thread Wolfgang Bangerth



Also another thing which gives me a headache is that why is the jacobian value
twice the value it should be. For a rectangular bilinear lagrange element we
can compute simply J by a/2 and b/2, but we would receive 1 and 1 for the
above case where a=b=2.
This doesn't correspond to what deal.II outputs.


It's not clear to me what you mean by a and b, and "1 and 1". Can you say what 
you expect the Jacobian matrix to be, and what deal.II outputs?




Additionally, using

std::vector > J = fe_values.get_jacobians ();

gives a vector containing n x n matrices?


Correct.


Since I can just output J[0][i][j], J[1][i][j] or J[10][i][j]. All work, but
why? Shouldn't it be somehow limited? (Stupid question maybe...)


What do you mean by "limited"? And do the J[0][i][j] match what you expect 
them to be?


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.


Re: [deal.II] Postprocessing on Vector Valued Problem

2016-12-21 Thread Wolfgang Bangerth

On 12/21/2016 02:05 PM, Jaekwang Kim wrote:


I tried to code following data post processing for my goal.


template

void

ShearProblem::post_processing ()

{



// std::cout << "   Post Processing..." << std::endl;



constMappingQ mapping (degree);



QGauss   quadrature_formula(1);



FEValues fe_values (mapping, fe, quadrature_formula,

 update_values|

 update_quadrature_points  |

 update_JxW_values |

 update_gradients);



constunsignedint  dofs_per_cell   = fe.dofs_per_cell;

constunsignedint  n_q_points  = quadrature_formula.size();



std::vector > local_symgrad_phi_u (n_q_points);

std::cout<< "   n_q_point: "<< n_q_points << std::endl;


constFEValuesExtractors::Vector velocities (0);



Vector cellwise_shear_rate (triangulation.n_active_cells());

// I want to save my data in this part and attach this to dof_handler,
so I can make data out in -vtk format later





unsignedintk=0;

typenameDoFHandler::active_cell_iterator

cell = dof_handler.begin_active(),

endc = dof_handler.end();

for(; cell!=endc; ++cell)

{

fe_values.reinit (cell);

fe_values[velocities].get_function_symmetric_gradients (solution,
local_symgrad_phi_u);





unsignedintq=0;



doubleshear_rate = std::sqrt( local_symgrad_phi_u[q]*
local_symgrad_phi_u[q] *2);


*cellwise_shear_rate(k)=shear_rate;*


k+=1;

}



}


I could able to - shear_rate at each cell appropriately,

 but How can I relate my 'cellwise_shear_rate' with dof handler so that I
finally able to see the values in 'vtk' format?
for now, it just seems that no more than a list of arrays


You can just attach this vector of values to a DataOut object and output that. 
That's what we often do with the vector of error indicators (which also has 
the size triangulation.n_active_cells()), but your case is really no different.


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] Postprocessing on Vector Valued Problem

2016-12-21 Thread Jaekwang Kim

Hi, all. 
I have a question on post processing of solution values. 
My solution is block vector solution that is consist of two parts. 
First part is velocities(u_x,u_y) and pressure. 

My purpose is to calculate shear-rate at each dof (or at each cell) which 
only deals with dof of velocities part. 

I tried to code following data post processing for my goal. 


template 

void

ShearProblem::post_processing ()

{



// std::cout << "   Post Processing..." << std::endl;



const MappingQ mapping (degree);



QGauss   quadrature_formula(1);



FEValues fe_values (mapping, fe, quadrature_formula,

 update_values|

 update_quadrature_points  |

 update_JxW_values |

 update_gradients);



const unsigned int   dofs_per_cell   = fe.dofs_per_cell;

const unsigned int   n_q_points  = quadrature_formula.size();



std::vector > local_symgrad_phi_u 
(n_q_points);

std::cout<< "   n_q_point: " << n_q_points << std::endl;


const FEValuesExtractors::Vector velocities (0);



Vector cellwise_shear_rate (triangulation.n_active_cells());

// I want to save my data in this part and attach this to 
dof_handler, so I can make data out in -vtk format later





unsigned int k=0;

typename DoFHandler::active_cell_iterator

cell = dof_handler.begin_active(),

endc = dof_handler.end();

for (; cell!=endc; ++cell)

{

fe_values.reinit (cell);

fe_values[velocities].get_function_symmetric_gradients 
(solution, local_symgrad_phi_u);





unsigned int q=0;



double shear_rate = std::sqrt( local_symgrad_phi_u[q]* 
local_symgrad_phi_u[q] *2 );


*cellwise_shear_rate(k)=shear_rate;*


k+=1;

}



}

I could able to - shear_rate at each cell appropriately,

 but How can I relate my 'cellwise_shear_rate' with dof handler so that I 
finally able to see the values in 'vtk' format? 
for now, it just seems that no more than a list of arrays


Always thank you all!

Jaekwang Kim 

-- 
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] Re: deal.II FE structure

2016-12-21 Thread Jean-Paul Pelteret
Dear Seyed,

You may also be interested in looking at the module on the interplay 
between finite elements and mappings 
, and the 
description to the mapping 
 class. Note 
that the reference element is [0,1]^{dim}.

Regards,
Jean-Paul

On Wednesday, December 21, 2016 at 8:32:31 PM UTC+1, Denis Davydov wrote:
>
> Hi Seyed,
>
> On Wednesday, December 21, 2016 at 7:16:32 PM UTC+1, mohseng...@gmail.com 
> wrote:
>>
>> Hi,
>>
>> I am trying to understand and learn the main structure of deal.II's 
>> design, because I cannot work with a program until I do not understand it 
>> fully.
>> Hence, I am checking a simple 2D cube with hand-written FE solutions to 
>> see how things are stored in deal.II.
>>
>> The model I use is created directly in deal.II by means of
>>
>> GridGenerator::hyper_cube (triangulation, -1, 1);
>> triangulation.refine_global (0);
>>
>> So a cube with a length of 2. 
>>
>> This is what I get for a 2D cube with linear shape functions and 4 GAUSS 
>> points (each row of the matrices represents a GP):
>>
>> GAUSS POINT COORDINATES
>>
>> -0.57735   -0.57735   
>> 0.57735   -0.57735   
>> -0.57735   0.57735   
>> 0.57735   0.57735   
>>
>>
>> SHAPE FUNCTIONS
>>
>>   0.6220   0.6220   0.1667   0.1667   0.1667   0.1667   0.0447   0.0447 
>>   0.1667   0.1667   0.6220   0.6220   0.0447   0.0447   0.1667   0.1667 
>>   0.1667   0.1667   0.0447   0.0447   0.6220   0.6220   0.1667   0.1667 
>>   0.0447   0.0447   0.1667   0.1667   0.1667   0.1667   0.6220   0.6220 
>>
>>
>> SHAPE FUNCTION DERIVATIVES
>>
>> -0.394338 -0.394338   -0.394338 -0.394338   0.394338 -0.105662   0.394338 
>> -0.105662   -0.105662 0.394338   -0.105662 0.394338   0.105662 0.105662   
>> 0.105662 0.105662   
>> -0.394338 -0.105662   -0.394338 -0.105662   0.394338 -0.394338   0.394338 
>> -0.394338   -0.105662 0.105662   -0.105662 0.105662   0.105662 0.394338   
>> 0.105662 0.394338   
>> -0.105662 -0.394338   -0.105662 -0.394338   0.105662 -0.105662   0.105662 
>> -0.105662   -0.394338 0.394338   -0.394338 0.394338   0.394338 0.105662   
>> 0.394338 0.105662   
>> -0.105662 -0.105662   -0.105662 -0.105662   0.105662 -0.394338   0.105662 
>> -0.394338   -0.394338 0.105662   -0.394338 0.105662   0.394338 0.394338   
>> 0.394338 0.394338   
>>
>>
>> JACOBIAN MATRIX
>>
>> 2 0 
>> 0 2 
>>
>> Unfortunately, I do not understand the point why the shape function 
>> derivatives change, if I create a geometry by means of:
>>
>> GridGenerator::hyper_cube (triangulation, 0, 1);
>> triangulation.refine_global (0);
>>
>
> It's because the gradient is w.r.t. the real space, as opposed to the 
> natural coordinates.
> Check out this lecture by Prof. Bangerth 
> http://www.math.colostate.edu/~bangerth/videos.676.10.html
> Actually I would recommend to go through ALL of them, and not just this 
> one, which answers this particular question.
>
> Regards,
> Denis.
>
>
>>

-- 
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] Re: deal.II FE structure

2016-12-21 Thread Denis Davydov
Hi Seyed,

On Wednesday, December 21, 2016 at 7:16:32 PM UTC+1, mohseng...@gmail.com 
wrote:
>
> Hi,
>
> I am trying to understand and learn the main structure of deal.II's 
> design, because I cannot work with a program until I do not understand it 
> fully.
> Hence, I am checking a simple 2D cube with hand-written FE solutions to 
> see how things are stored in deal.II.
>
> The model I use is created directly in deal.II by means of
>
> GridGenerator::hyper_cube (triangulation, -1, 1);
> triangulation.refine_global (0);
>
> So a cube with a length of 2. 
>
> This is what I get for a 2D cube with linear shape functions and 4 GAUSS 
> points (each row of the matrices represents a GP):
>
> GAUSS POINT COORDINATES
>
> -0.57735   -0.57735   
> 0.57735   -0.57735   
> -0.57735   0.57735   
> 0.57735   0.57735   
>
>
> SHAPE FUNCTIONS
>
>   0.6220   0.6220   0.1667   0.1667   0.1667   0.1667   0.0447   0.0447 
>   0.1667   0.1667   0.6220   0.6220   0.0447   0.0447   0.1667   0.1667 
>   0.1667   0.1667   0.0447   0.0447   0.6220   0.6220   0.1667   0.1667 
>   0.0447   0.0447   0.1667   0.1667   0.1667   0.1667   0.6220   0.6220 
>
>
> SHAPE FUNCTION DERIVATIVES
>
> -0.394338 -0.394338   -0.394338 -0.394338   0.394338 -0.105662   0.394338 
> -0.105662   -0.105662 0.394338   -0.105662 0.394338   0.105662 0.105662   
> 0.105662 0.105662   
> -0.394338 -0.105662   -0.394338 -0.105662   0.394338 -0.394338   0.394338 
> -0.394338   -0.105662 0.105662   -0.105662 0.105662   0.105662 0.394338   
> 0.105662 0.394338   
> -0.105662 -0.394338   -0.105662 -0.394338   0.105662 -0.105662   0.105662 
> -0.105662   -0.394338 0.394338   -0.394338 0.394338   0.394338 0.105662   
> 0.394338 0.105662   
> -0.105662 -0.105662   -0.105662 -0.105662   0.105662 -0.394338   0.105662 
> -0.394338   -0.394338 0.105662   -0.394338 0.105662   0.394338 0.394338   
> 0.394338 0.394338   
>
>
> JACOBIAN MATRIX
>
> 2 0 
> 0 2 
>
> Unfortunately, I do not understand the point why the shape function 
> derivatives change, if I create a geometry by means of:
>
> GridGenerator::hyper_cube (triangulation, 0, 1);
> triangulation.refine_global (0);
>

It's because the gradient is w.r.t. the real space, as opposed to the 
natural coordinates.
Check out this lecture by Prof. 
Bangerth http://www.math.colostate.edu/~bangerth/videos.676.10.html
Actually I would recommend to go through ALL of them, and not just this 
one, which answers this particular question.

Regards,
Denis.


>

-- 
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] deal.II FE structure

2016-12-21 Thread mohsengineering
Hi,

I am trying to understand and learn the main structure of deal.II's design, 
because I cannot work with a program until I do not understand it fully.
Hence, I am checking a simple 2D cube with hand-written FE solutions to see 
how things are stored in deal.II.

The model I use is created directly in deal.II by means of

GridGenerator::hyper_cube (triangulation, -1, 1);
triangulation.refine_global (0);

So a cube with a length of 2. 

This is what I get for a 2D cube with linear shape functions and 4 GAUSS 
points (each row of the matrices represents a GP):

GAUSS POINT COORDINATES

-0.57735   -0.57735   
0.57735   -0.57735   
-0.57735   0.57735   
0.57735   0.57735   


SHAPE FUNCTIONS

  0.6220   0.6220   0.1667   0.1667   0.1667   0.1667   0.0447   0.0447 
  0.1667   0.1667   0.6220   0.6220   0.0447   0.0447   0.1667   0.1667 
  0.1667   0.1667   0.0447   0.0447   0.6220   0.6220   0.1667   0.1667 
  0.0447   0.0447   0.1667   0.1667   0.1667   0.1667   0.6220   0.6220 


SHAPE FUNCTION DERIVATIVES

-0.394338 -0.394338   -0.394338 -0.394338   0.394338 -0.105662   0.394338 
-0.105662   -0.105662 0.394338   -0.105662 0.394338   0.105662 0.105662   
0.105662 0.105662   
-0.394338 -0.105662   -0.394338 -0.105662   0.394338 -0.394338   0.394338 
-0.394338   -0.105662 0.105662   -0.105662 0.105662   0.105662 0.394338   
0.105662 0.394338   
-0.105662 -0.394338   -0.105662 -0.394338   0.105662 -0.105662   0.105662 
-0.105662   -0.394338 0.394338   -0.394338 0.394338   0.394338 0.105662   
0.394338 0.105662   
-0.105662 -0.105662   -0.105662 -0.105662   0.105662 -0.394338   0.105662 
-0.394338   -0.394338 0.105662   -0.394338 0.105662   0.394338 0.394338   
0.394338 0.394338   


JACOBIAN MATRIX

2 0 
0 2 

Unfortunately, I do not understand the point why the shape function 
derivatives change, if I create a geometry by means of:

GridGenerator::hyper_cube (triangulation, 0, 1);
triangulation.refine_global (0);

Does this mean that we define the cube in natural coordinates? From what I 
know, the shape functions only depend on natural coordinates which means 
independent from the geometry, so they should remain the same, if the 
natural coordinates from -1 to 1 remain.
Also another thing which gives me a headache is that why is the jacobian 
value twice the value it should be. For a rectangular bilinear lagrange 
element we can compute simply J by a/2 and b/2, but we would receive 1 and 
1 for the above case where a=b=2.
This doesn't correspond to what deal.II outputs.

Additionally, using 

std::vector > J = fe_values.get_jacobians ();
 
gives a vector containing n x n matrices?

Since I can just output J[0][i][j], J[1][i][j] or J[10][i][j]. All work, 
but why? Shouldn't it be somehow limited? (Stupid question maybe...)


I assume such simple questions may annoy the deal.II developers, I know. 
But I promise whenever I learn the program thoroughly I come up with more 
complex questions ;)


Best regards,
Seyed Ali 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.


[deal.II] restart/checkpoint strategies in user applications/libraries

2016-12-21 Thread Denis Davydov
Hi all,

This is a bit off-topic, but i would like to ask how do you prefer to 
implement restart/checkpoints in your user applications?
I see two options:

1. Make restarting usable with the same input parameter file so that 
calculations will continue exactly from the point of restart.
For example, if you wanted to do 10 adaptive refinements and one node 
failed on 6th step, you would use exactly the same input parameter 
file to automatically continue from the 5th step if the application see 
that restart information is available.

2. Another approach is to consider restarts as an arbitrary input (similar 
to the input mesh / initial conditions and alike).
With the same example in mind, if one uses the same input parameter file 
with 10 refinement steps
but specify `Restart = true` and thereby use the restart generated from 5th 
step,
then the program would actually end up doing 10 refinement steps starting 
from the input data taken from the restart.
In other words, 15 steps in total.

Maybe ASPECT (Timo, Wolfagng?), pi-DoMUS ( Luca?) or DOpElib developers 
could comment on which restart strategy you have chosen 
for your libraries and why?

Regards,
Denis.

-- 
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.