Ernesto Ismail
I am using deal.II for a problem where the variable of interest is a second order Symmetric Tensor (e.g solving for the stress in an evolution equation with a given velocity field). My question is in two parts, both relating to derivatives of Symmetric Tensors.
1) I have a field of symmetric tensors that I have found. I need to use the divergence of this field in a second computation. Were it a vector field I would use get_function_gradients and then take the trace of this for the divergence. How should I go about doing this for the SymmetricTensor<2,dim> field?
2) I am building an equation where I am solving for a field of
symmetric tensors. At assembly I would like to use the gradient
associated with the shape functions of the elements describing them. For
a vector valued problem there is a "view" type which allows me to
extract the symmetric tensor itself (const Tensor<2,dim> phi_j_T
= fe_values[Ts].value (j, q);) but there doesn't seem to be an
equivalent (const Tensor<3,dim> phi_j_T_grad =
fe_values[Ts].gradient (j, q);) How should I tackle this?
Thanks in advance,
Ernesto
Wolfgang Bangerth
> 1) I have a field of symmetric tensors that I have found. I need to use the
> divergence of this field in a second computation. Were it a vector field I
> would use get_function_gradients and then take the trace of this for the
> divergence. How should I go about doing this for the SymmetricTensor<2,dim> field?
FEValuesExtractors::SymmetricTensor<2> Ts;
fe_values[Ts].get_function_divergences (input, output);
> 2) I am building an equation where I am solving for a field of symmetric
> tensors. At assembly I would like to use the gradient associated with the
> shape functions of the elements describing them. For a vector valued problem
> there is a "view" type which allows me to extract the symmetric tensor itself
> (const Tensor<2,dim> phi_j_T = fe_values[Ts].value (j, q);) but there
> doesn't seem to be an equivalent (const Tensor<3,dim> phi_j_T_grad =
> fe_values[Ts].gradient (j, q);) How should I tackle this?
What code do you have that doesn't compile for you?
Best
W.
--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@math.tamu.edu
www: http://www.math.tamu.edu/~bangerth/
Ernesto Ismail
With regards to point 1)
I had not noted that there was a get_function_divergences function:
https://www.dealii.org/developer/doxygen/deal.II/classFEValuesViews_1_1Vector.html#aa20eae586e7a296f95770cb02c5ad773
However, a function to find the gradients does not seem to exist eg if I use:
fe_values[Ts].get_function_gradients (old_solution, old_T_grads);
I get the compiler error:
error: ‘const class dealii::FEValuesViews::SymmetricTensor<2, 2, 2>’ has no member named ‘get_function_gradients’
With regards to point 2)
I am trying to do something along the lines of:
const Tensor<3,dim> phi_i_P_grad = fe_values[Ps].gradient (i, q);
The compiler error I get is:
error: ‘const class dealii::FEValuesViews::SymmetricTensor<2, 2, 2>’ has no member named ‘gradient’
Am I correct in my understanding that at present it is not possible to find gradients of second order tensors directly? If yes, how would I go about finding them?
Regards,
Ernesto
Wolfgang Bangerth
> However, a function to find the gradients does not seem to exist eg if I use:
> fe_values[Ts].get_function_gradients (old_solution, old_T_grads);
>
> I get the compiler error:
> error: ‘const class dealii::FEValuesViews::SymmetricTensor<2, 2, 2>’ has no
> member named ‘get_function_gradients’
>
> With regards to point 2)
> I am trying to do something along the lines of:
> const Tensor<3,dim> phi_i_P_grad = fe_values[Ps].gradient (i, q);
>
> The compiler error I get is:
> error: ‘const class dealii::FEValuesViews::SymmetricTensor<2, 2, 2>’ has no
> member named ‘gradient’
>
> Am I correct in my understanding that at present it is not possible to find
> gradients of second order tensors directly? If yes, how would I go about
> finding them?
get_function_gradients() functions in this class.
I think the simplest way for you around this problem would actually be to add
these missing functions to the FEValuesViews::SymmetricTensor class. Take a
look at include/deal.II/fe/fe_values.h, starting around line 900 of the
current development sources. The function already has divergence() and
get_function_divergences() members. It shouldn't be too difficult to write the
missing functions in a similar style. The computed type should in both cases
be Tensor<3,spacedim>.
I think adding these functions is going to be simpler than trying to work
around the lack thereof.
Of course we would be very happy to take whatever code you write. Let us know
if you need help!
Andrew Mcbride
You might want to make the return type a SymmetricTensor<3,spacedim> as you will have symmetry in A_ij,k = A_ji,k as A_ij = A_ji. That should be fine as long as you define symmetry of a third-order tensor to be in the first two slots.
Wolfgang is correct, the functions are not there. The reason for this was that we simply did not need them at the time.
Andrew
> 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+un...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
Andrew Mcbride
maybe stick with Wolfgang’s suggestion
A
Ernesto Ismail
I will take a look at the existing FEValuesViews::SymmetricTensor class and try write the functions. I'll try get to this early next week.
You received this message because you are subscribed to a topic in the Google Groups "deal.II User Group" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/dealii/79newCZGhhM/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dealii+un...@googlegroups.com.
Wolfgang Bangerth
>
> From my understanding of the SymmetricTensor class it only supports
> orders 2,4 etc. The wording from the documentation is "Provide a class
> that stores symmetric tensors of rank 2,4,..." from
> https://www.dealii.org/developer/doxygen/deal.II/classSymmetricTensor.html
> I will take a look at the existing FEValuesViews::SymmetricTensor class
> and try write the functions. I'll try get to this early next week.
Make sure you start from the current version you can find on github as
that makes later merging the changes a lot easier!
Best
Wolfgang
Ernesto Ismail
Ernesto
--
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
Wolfgang Bangerth
> I just pulled the latest version and am planning my approach. Once I'm
> done, how would you like me to contribute it? Is it possible for me to
> push my submission to the repository for review?
name; you then pull a copy of this fork onto your local harddrive (i.e.,
you "clone" it). You should then make all of your changes on a branch,
push this branch to your forked repository, and then create a "pull
request" that we can look at, comment on, and ultimately merge.
There are a number of tutorials online that illustrate how this works if
you are unfamiliar with this language.
Best
W.
Metin C.
Wolfgang Bangerth
Metin,
> I need help to figure out some functions in fe_values.h and fe_values.cc
> sources.
>
> Could you please help me understand what's happening in the divergence
> function (SymmetricTensor<2, dim, spacedim>::divergence) and the
> symmetric_gradient function for the vectors
> (Vector<dim,spacedim>::symmetric_gradient), or point out any material
> (manual, video etc.)?
have trouble with?
Best
Wolfgang
--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/
Metin Cakircali
for (unsigned int shape_function=0;
shape_function<dofs_per_cell; ++shape_function)
{
for (unsigned int d=0; d<spacedim; ++d)
if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
{
const dealii::Tensor<1,spacedim> *shape_gradient_ptr =
&shape_gradients[shape_function_data[shape_function].
row_index[d]][0];
for (unsigned int q_point=0; q_point<n_quadrature_points; ++q_point)
divergences[q_point] += value * (*shape_gradient_ptr++)[d];
}
}
To me, it's not clear how the "value and shape_gradient_ptr" is used or in general how the (e.g.) u_x+v_y+w_z is done..
For the SymmetricTensor, the code differs as;
const unsigned int ii = dealii::SymmetricTensor<2,spacedim>::
unrolled_to_component_indices(comp)[0];
const unsigned int jj = dealii::SymmetricTensor<2,spacedim>::
unrolled_to_component_indices(comp)[1];
for (unsigned int q_point = 0; q_point < n_quadrature_points;
++q_point, ++shape_gradient_ptr)
{
divergences[q_point][ii] += value * (*shape_gradient_ptr)[jj];
if (ii != jj)
divergences[q_point][jj] += value * (*shape_gradient_ptr)[ii];
}
I again cannot follow it. and what should I understand, that the loop over spacedim (or SymmetricTensor<2,spacedim>::n_independent_components) not implemented or not needed? please see, Assert (false, ExcNotImplemented());
Also, as a general question, how/where does the "Tensor<2, dim, spacedim>::divergence" function operate?
Daniel Arndt
[...]
For the vectors, the implementation (for the else case) is as following;for (unsigned int shape_function=0;
shape_function<dofs_per_cell; ++shape_function)
{
for (unsigned int d=0; d<spacedim; ++d)
if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
{
const dealii::Tensor<1,spacedim> *shape_gradient_ptr =
&shape_gradients[shape_function_data[shape_function].
row_index[d]][0];
for (unsigned int q_point=0; q_point<n_quadrature_points; ++q_point)
divergences[q_point] += value * (*shape_gradient_ptr++)[d];
}
}
To me, it's not clear how the "value and shape_gradient_ptr" is used or in general how the (e.g.) u_x+v_y+w_z is done..
Metin Cakircali
if (snc != -1)
{
const unsigned int comp =
shape_function_data[shape_function].single_nonzero_component_index;
const dealii::Tensor < 1, spacedim> *shape_gradient_ptr =
&shape_gradients[snc][0];
const TableIndices<2> indices = dealii::Tensor<2,spacedim>::unrolled_to_component_indices(comp);
const unsigned int ii = indices[0];
const unsigned int jj = indices[1];
for (unsigned int q_point = 0; q_point < n_quadrature_points;
++q_point, ++shape_gradient_ptr)
{
divergences[q_point][jj] += value * (*shape_gradient_ptr)[ii];
}
}
else
{
for (unsigned int d = 0;
d < dim*dim; ++d)
if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
{
Assert (false, ExcNotImplemented());
}
}
There is only the "snc != 1" case, and it seems that the component now means the tensor components and not the directions in space. Is this function (for tensors) not yet implemented? or is it sth related to non-primitive shape functions (I would appreciate if you could comment on this as well)?
Thank you in advance.
Best,
Metin
Daniel Arndt
I understand that the components are in space (x/y/z), and (*shape_gradient_ptr++)[d] would be derivative of phi in each direction; is that correct?
(compared to vectors) the implementation of divergence function for the tensors is as following;if (snc != -1)
{
const unsigned int comp =
shape_function_data[shape_function].single_nonzero_component_index;
const dealii::Tensor < 1, spacedim> *shape_gradient_ptr =
&shape_gradients[snc][0];
const TableIndices<2> indices = dealii::Tensor<2,spacedim>::unrolled_to_component_indices(comp);
const unsigned int ii = indices[0];
const unsigned int jj = indices[1];
for (unsigned int q_point = 0; q_point < n_quadrature_points;
++q_point, ++shape_gradient_ptr)
{
divergences[q_point][jj] += value * (*shape_gradient_ptr)[ii];
}
}
There is only the "snc != 1" case, and it seems that the component now means the tensor components and not the directions in space. Is this function (for tensors) not yet implemented? or is it sth related to non-primitive shape functions (I would appreciate if you could comment on this as well)?