How to calculate the mean curvature

[Oded Yaakobi]

Dear all,

I would like to calculate the mean curvature of the boundary of a three-dimensional domain.
Does anybody know if and how it could be done in deal.ii?

Thanks,
Oded

[Wolfgang Bangerth]

On 05/29/2015 04:27 PM, Oded Yaakobi wrote:
>
> I would like to calculate the mean curvature of the boundary of a
> three-dimensional domain.
> Does anybody know if and how it could be done in deal.ii?

Can you state what the formula is for the thing you want to compute? This
would make it easier to suggest how to do it.

Best
W.

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

[Oded Yaakobi]

Dear Wolfgang,

 

There are various equivalent formulas for calculating the mean curvature. Some of them appear in the attached paper where it is denoted by K_M:

 

1) Page 635 - at the bottom of the page.

2) Page 642 - Eq. 4.2.

3) Page 645 - Corollary 4.5.

 

Thank you for your help,

Oded

[Wolfgang Bangerth]

Oded,

> There are various equivalent formulas for calculating the mean curvature. Some
> of them appear in the attached paper where it is denoted by K_M:
>
> 1) Page 635 - at the bottom of the page.
>
> 2) Page 642 - Eq. 4.2.
>
> 3) Page 645 - Corollary 4.5.

All of these formulas assume that your surface is smooth, but the surface of a
finite element mesh is only piecewise polynomial and, consequently, none of
these representations of the mean curvature apply. You probably want to figure
out how to suitably generalize the concept before looking for a way to
implement it.

[Oded Yaakobi]

Hi Wolfgang,

Thanks, I understand the issue that you have noted.

Instead of calculating the mean curvature, I can choose another way to tackle my problem. It requires calculating the gradient of the surface area of each face that is part of the boundary. This quantity should be calculated at the quadrature points of the boundary surface element, and the gradient should be taken with respect to the real coordinates of the corresponding quadrature point. Multiplying this quantity by a surface tension coefficient gives a load on each quadrature point. Then, a total surface tension force which is related to the face has to be calculated by taking into account the weight of each surface quadrature point.    

Is it possible to perform this calculation using deal.ii? If yes, could you please advise me how to do it?
In particular, I do not know how to calculate the aforementioned surface area gradient with respect to the real coordinates of the quadrature points.

Best,
Oded

[Wolfgang Bangerth]

Oded,

> Instead of calculating the mean curvature, I can choose another way to tackle
> my problem. It requires calculating the gradient of the surface area of each
> face that is part of the boundary. This quantity should be calculated at the
> quadrature points of the boundary surface element, and the gradient should be
> taken with respect to the real coordinates of the corresponding quadrature
> point. Multiplying this quantity by a surface tension coefficient gives a load
> on each quadrature point. Then, a total surface tension force which is related
> to the face has to be calculated by taking into account the weight of each
> surface quadrature point.

This is again a case where formulas would help :-) The description is too
complicated for me to follow.

I think whatever approach you come up with, you should think about how it
would apply to a domain that is described by linear triangles. For these, the
boundary consists of line segments with no curvature, and the traditional
definition of the curvature would yield delta functions at the vertices.

[Oded Yaakobi]

Hi Wolfgang,

I tried to be as concise as possible with my postings but obviously I had to elaborate more in order to clarify my intention.  

Attached is a more detailed description including formulas of my post from yesterday. If this approach is feasible, then there is no need to calculate the curvature.

Thank you very much for your help!
Oded

[Wolfgang Bangerth]

Oded,

> I tried to be as concise as possible with my postings but obviously I had to
> elaborate more in order to clarify my intention.
>
> Attached is a more detailed description including formulas of my post from
> yesterday. If this approach is feasible, then there is no need to calculate
> the curvature.

The term dA/dr is of course somehow related to the mapping from the reference
cell to the real cell. You'd have to think a bit about how they are related,
but it shouldn't be very hard.

There is a good amount of literature on terms like the ones you want to
consider. What do others do? I would suggest to also check out the papers by
my colleague Andrea Bonito who has several applications where surface tension
matters.

[Andrew Mcbride]

You might like to have a look at http://www.cerecam.uct.ac.za/code/surface_energy/doc/html/ and arXiv:1506.01361 

It includes an implementation of the slightly exotic case of surface tension for solids. You will need to read some of the references to see how the curvature is captured in the surface divergence operator.

A

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

[Oded Yaakobi]

Hi Wolfgang,

I figured out how could I calculate the required term dA/dr. I am working now on implementing it in DEAL-II. In order to do so, I have to calculate the shape function gradients of the transformation of the boundary face elements. I tested my code by modifying step-20 and adding the flag "update_transformation_gradients" to the variable "fe_face_values". However, when I try to read the value of this expression, the compiler gives an error:

error: ‘class dealii::FEFaceValues<3, 3>’ has no member named ‘transformation_gradients’

I guess that I don't use the correct syntax, but I weren't able to find out what the correct syntax should be. Could you please advise me how to extract the required values?

Attached is my slightly modified version of step-20. The error is related to
line 525: fe_face_values.transformation_gradients;

In addition, what should be the variable type of dSi_dxi below?
line 526: dSi_dxi = fe_face_values.transformation_gradients;

Thank you,
Oded

[Sungho Yoon]

Well. I read pdf files. it's still not clear. Which would you  use the implicit function or surface on the finite element mesh explicitly? If you are with implicit function with interface, worthy to find "stephane zaleski". He's a singular person about this.

[Oded Yaakobi]

Hi Sungho,

Thank you for your comment. I will take a look at Zaleski's publications, which are somewhat related to my research.

However, the issue that I was referring to in my post yesterday is not a mathematical difficulty. It is purely a question about how to use DEAL-II. As I noted, I tried to compile a slightly modified version of step-20 (which I attached to my post), but encountered a problem in doing so.

Best,
Oded

[Wolfgang Bangerth]

On 06/22/2015 03:09 PM, Oded Yaakobi wrote:
> Hi Wolfgang,
>
> I figured out how could I calculate the required term dA/dr. I am working now
> on implementing it in DEAL-II. In order to do so, I have to calculate the
> shape function gradients of the transformation of the boundary face elements.
> I tested my code by modifying step-20 and adding the flag
> "update_transformation_gradients" to the variable "fe_face_values". However,
> when I try to read the value of this expression, the compiler gives an error:
>
> error: ‘class dealii::FEFaceValues<3, 3>’ has no member named
> ‘transformation_gradients’

Correct. I needed to spend a bit looking through our codes to understand what
is actually happening, but the way I read the code is that this flag is (the
way it currently is) intended only for internal use: finite elements use it to
signal to the mappings that they need the derivatives of the mapping, for
example to compute co- or contravariant mappings. There is currently no
provision to copy this data from the Mapping internal data structures to
FEValues internal data structures, and there is also no function to allow
access to this information from the FEValues object.

Neither would be terribly difficult to implement. Let me know if you want to
go down that route and I can give you a few pointers to what needs to be changed.

[Oded Yaakobi]

Hi Wolfgang,

Thanks for clarifying this point. Yes, please let me know what should I do in order to implement the functionality that is currently missing in DEAL-II.

Best,

Oded

[Wolfgang Bangerth]

On 07/01/2015 09:50 AM, Oded Yaakobi wrote:
>
> Thanks for clarifying this point. Yes, please let me know what should I do in
> order to implement the functionality that is currently missing in DEAL-II.

I think you should try to understand how the corresponding update_jacobians
flag is processed, and how this information is propagated. In both cases, the
data is computed by the mappings. In your case, this data is only used
internally, whereas in the case of the Jacobians, there must be a place where
the mapping copies the Jacobians into a structure provided by FEValues, and
then there are FEValues functions that allow the user to access the data.

You need to replicate these steps, i.e., the FEValues object must provide a
place similar to where it stores the Jacobians, the mappings must be able to
copy their internal data into these places, and then FEValues must gain
functions that allow to access this data.

Let me know if you have more concrete questions.

[Oded Yaakobi]

Hi Wolfgang,

I read more thoroughly the library documentation and I am little bit confused now. I am interested to get the derivatives of the shape functions of a face of a 3D volume element with respect to the natural coordinates of the reference 2D cell. Isn't that exactly the output of one of the functions or one of the attributes that are listed below?
If so, can't I use it directly instead of modifying the DEAL-II library source files?

By the way, I don't know whether it makes any difference, but my actual project is a variation of the Stokes problem. Hence I am using step-22 as a base to start with, where the only type of finite elements is FE_Q. 

Best,
Oded


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

template<int dim, int spacedim = dim>
virtual Tensor<1,dim> FiniteElement< dim, spacedim >::shape_grad     (     const unsigned int      i,
        const Point< dim > &      p
    )         const

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

template<int dim, int spacedim = dim>
virtual Tensor<1,dim> FiniteElement< dim, spacedim >::shape_grad_component     (     const unsigned int      i,
        const Point< dim > &      p,
        const unsigned int      component
    )         const
    virtual

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

template<class POLY, int dim = POLY::dimension, int spacedim = dim>
std::vector<std::vector<Tensor<1,dim> > > FE_Poly< POLY, dim, spacedim >::InternalData::shape_gradients

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

template<int dim, int spacedim = dim>
std::vector<Tensor<1,dim> > MappingQ1< dim, spacedim >::InternalData::shape_derivatives

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

template<int dim, int spacedim = dim>
virtual Tensor<1,dim> FESystem< dim, spacedim >::shape_grad     (     const unsigned int      i,
        const Point< dim > &      p
    )         const
    virtual

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[Wolfgang Bangerth]

Oded,

> I read more thoroughly the library documentation and I am little bit confused
> now. I am interested to get the derivatives of the shape functions of a face
> of a 3D volume element with respect to the natural coordinates of the
> reference 2D cell. Isn't that exactly the output of one of the functions or
> one of the attributes that are listed below?

Yes. But your previous question was about the derivative of the mapping from
reference to real cell, or did I completely misunderstand?

[Oded Yaakobi]

Wolfgang,

My question from last week aimed to get the same quantity, i.e. the derivatives of the shape functions of a face of a 3D volume element with respect to the natural coordinates of the reference 2D cell.

I thought that this is what would be obtained from the class "FEFaceValues" by turning on the flag "update_transformation_gradients". In the documentation of "Finite element access/FEValues classes" the description of this flag is "Shape function gradients of transformation. Compute the shape function gradients of the transformation defined by the Mapping." 

See:

http://dealii.org/developer/doxygen/deal.II/group__feaccess.html#ggaa94b67d2fdcc390690c523f28019e52fa78349f26f4f9386c740940ddbc4be585

I thought that this would give the derivatives of the shape functions with respect to the natural coordinates. This is in contrast to the flag "update_gradients" whose description at the same place is "Shape function gradients. Compute the gradients of the shape functions in coordinates of the real cell."

However, now it seems to me that I didn't understand properly the description of the flag "update_transformation_gradients", since the quantity the I am interested in has nothing to do with the mapping (although I still don't understand what is the meaning of this flag).

Anyway, if I understand correctly your last reply, then one (or more) of the functions that I listed yesterday is the one that I need. Which one of them should I use? 

If for instance I would like to add a call for this function into step-22, what would be the exact syntax of this command? 

Thank you very much for your continuous help,

Oded

[Wolfgang Bangerth]

Oded,

> My question from last week aimed to get the same quantity, i.e. the
> derivatives of the shape functions of a face of a 3D volume element with
> respect to the natural coordinates of the reference 2D cell.

Ah, then I misunderstood. So this is information that only involves the finite
element (as it does not involve any mapping) and that you can get directly
from the FE object. No FEValues object is necessary for that.

> I thought that this is what would be obtained from the class "FEFaceValues" by
> turning on the flag "update_transformation_gradients". In the documentation of
> "Finite element access/FEValues classes" the description of this flag is
> "Shape function gradients of transformation. Compute the shape function
> gradients of the transformation defined by the Mapping."
>
> See:
>
> http://dealii.org/developer/doxygen/deal.II/group__feaccess.html#ggaa94b67d2fdcc390690c523f28019e52fa78349f26f4f9386c740940ddbc4be585
>
> I thought that this would give the derivatives of the shape functions with
> respect to the natural coordinates. This is in contrast to the flag
> "update_gradients" whose description at the same place is "Shape function
> gradients. Compute the gradients of the shape functions in coordinates of the
> real cell."

The key word in the description of the former is *mapping*. It refers to the
mapping from the reference to the real cell. You don't seem to be interested
in that. But I think you already found that out.

> Anyway, if I understand correctly your last reply, then one (or more) of the
> functions that I listed yesterday is the one that I need. Which one of them
> should I use?

FiniteElement::shape_gradient for scalar elements. shape_gradient_component is
for vector-valued elements. The finite element you operate on is, I suppose,
the one that you use for your discretization.