Defining tensor in different dim and spacedim

[benhour....@gmail.com]

Dear All,

I am solving an elastic problem coupled with heat equation transfer (thermal-elasticity). I should solve the problem in a curved domain (circle in 2d and sphere in 3d) . I am developing step 44 in this case. I have to define some symmetric tensors. I have attached some parts of my code that i have problem with. I have two unknown variables. one for displacement (vector) and another for other parameter(scalar). I do really appreciate your kindness if you let me know how I can define symmetric tensor or in general tensor for this case. In addition, should I define DofHandler<dim, spacedim> for both scalar and vector valued unknown? My main problem is how to define tensor in different dim and spacedim problem.

       template <int spacedim>

       class PointHistory

          {

           public:

           PointHistory()

             :

           material(NULL),

      F_inv(StandardTensors<spacedim>::I),

              elasStress(SymmetricTensor<2, spacedim>()),

      totalStress(SymmetricTensor<2, spacedim>()),

              Jc(SymmetricTensor<4, spacedim>())

         {}

        virtual ~PointHistory()

          {

            delete material;

            material = NULL;

          }

        void setup_lqp (const Parameters::AllParameters &parameters)

          {

            material = new Material_Variables<spacedim>(parameters.lambda0,

            parameters.mu0, parameters.lambda1, parameters.mu1, parameters.chi, parameters.H,

            parameters.thetain, parameters.thetacr, parameters.thetaeq, parameters.theta, parameters.gammasv);

            update_values(Tensor<2, dim>(), double ());

          }

        void update_values (const Tensor<2, dim> &Grad_u_n, const double eta

            )

          {

            const Tensor<2, dim> eps = 0.5 * (Grad_u_n + transpose(Grad_u_n));                                            // Total Strain

            const double phi = std::pow(eta, 2) * (3 - 2 * eta);                                                                         // Function of order parameter

            dphi = 6 * eta * (1 - eta);                                                                                                             // Derivative of the interpolation function

            const double doubleWellBarrier = 2*eta + 4 * std::pow(eta, 3) - 6*std::pow(eta, 2);

            const double etaBVG = -6 * chi1 * delg1 * (1+eps_V) * dphi;                                                         // Boundary value for G-L equation

            const double etaBVM = 2 * gam1 / 30e-9; // Boundary value for Mechanics

// Thermal Strain Tensor

       const SymmetricTensor<2, dim> epsThermal = alphas * (thetaeq1 - thetain1) * StandardTensors<dim>::I +

       (alpham + (alphas - alpham) * phi) * (theta1 - thetaeq1) * StandardTensors<dim>::I;

// Transformation Strain Tensor

       const SymmetricTensor<2, dim> epsTransformation = (1.0/3.0) * 0.02 * StandardTensors<dim>::I *

                                               (1 - phi);

//The Initial Strain

       const SymmetricTensor<2, dim> epsInitial = SymmetricTensor<2, dim>(epsThermal) +

        (1-phi) * 0.02 * StandardTensors<dim>::I;

// Elastic Strain Tensor

       const SymmetricTensor<2, dim> epsElastic = symmetrize(Tensor<2, dim>(eps)-

                                          SymmetricTensor<2, dim>(epsTransformation)-

                                          SymmetricTensor<2, dim>(epsInitial) -

          SymmetricTensor<2, dim>(epsThermal));

// The Deviatoric Strain

const SymmetricTensor<2, dim> epsDeviatoric = symmetrize(Tensor<2, dim>(eps)-(1.0/3.0) * eps_V*

           StandardTensors<dim>::I);

       material->update_material_data(F, eps, phi, dphi, doubleWellBarrier, etaBVG, etaBVM

        );

       F_inv = invert(F);

       eps_V = material->get_eps_V();

       elasStress = material->get_elasStress();

//        surfTenStress = material->get_surfTenStress();   // It is exactly equal to initial stress...

       totalStress = material->get_totalStress();       // Extracting total Cauchy stress

       Jc = material->get_Jc();                                             // Extracting elasticity stiffnes tensor

       driving_force_noStress = material->get_driving_force_noStress();   // Extracting driving force with no stress

       eta_boundary_valueG = material->get_eta_boundary_valueG();                                      // Extracting the boundary value for G-L equation

                eta_boundary_valueM = material->get_eta_boundary_valueM();                                       // Extracting the boundary value for mechanics

       const double temp = -(1.0/(1+eps_V))*0.02*mean_elas_stress;

       const double temp1 = temp + 3*(1.0/(1+eps_V))*mean_elas_stress*

                               (alphas - alpham)*(theta1 - thetaeq1);

       const double temp2 = temp1 - ((1.0/(1+eps_V))*(0.5*(lambda1-lambda0)*

                              SymmetricTensor<2, dim>(epsElastic) * SymmetricTensor<2, dim>(epsElastic)));

       const double temp3 = -(1.0/(1+eps_V))*mu*SymmetricTensor<2, dim>(epsDeviatoric)*

                               SymmetricTensor<2, dim>(epsDeviatoric);

// Computing total driving force

       get_driv_forc = temp2 * temp3 * dphi + driving_force_noStress;

// computing the total boundary value for the G-L equation

            get_eta_bou_valG = eta_boundary_valueG;

// computing the total boundary value for the Mechanics

            get_eta_bou_valM = eta_boundary_valueM;

       Assert(determinant(F_inv) > 0, ExcInternalError());

        }

        double get_det_F() const

        {

          return material->get_det_F();

        }

        double get_eps_V() const

        {

             return material->get_eps_V();

        }

       const Tensor<2, dim> &get_F_inv() const

       {

          return F_inv;

       }

        const SymmetricTensor<2, dim> &get_elasStress() const

{

        return elasStress;

}

        const SymmetricTensor<2, dim> &get_totalStress() const

        {

            return totalStress;

        }

//        const SymmetricTensor<2, dim> &get_surfTenStress() const

// {

//         return surfTenStress;

// }

        const SymmetricTensor<4, dim> &get_Jc() const

        {

            return Jc;

        }

        double get_driving_force() const

        {

            return get_driv_forc;

        }

        double get_meanElStress() const

        {

        return mean_elas_stress;

        }

// The boundary value for the G-L equation

        double get_eta_orderG() const

        {

        return get_eta_bou_valG;

        }

// The boundary value for the Mechanics

        double get_eta_orderM() const

{

        return get_eta_bou_valM;

}

        private:

        Material_Variables<dim>             *material;

        Tensor<2, dim>                           eps;

        Tensor<2, dim>                           F_inv;

        SymmetricTensor<2, dim>           epsTransformation;

        SymmetricTensor<2, dim>           epsThermal;

        SymmetricTensor<2, dim>           epsInitial;

        SymmetricTensor<2, dim>           epsElastic;

        SymmetricTensor<2, dim>           epsDeviatoric;

        double                                        eps_V;

        Tensor<2, dim>                           elasStress;

        Tensor<2, dim>                           totalStress;

//      SymmetricTensor<2, dim>           initialStress;

        SymmetricTensor<4, dim>           Jc;

[Jean-Paul Pelteret]

Dear Benhour,

I'd like to help you but I really don't understand your question. Since its templated on the [space] dimension, the PointHistory class that you've displayed here is potentially suitable for either codimension 0 problems (dependent on your implementation of a constitutive law, of course, and assuming that you swap all instances of "dim" for the template parameter "spacedim"). Can you please try to explain a little more clearly what you are trying to do and what the issue is, and give some more context to the code you've quoted?

- Are you solving your problem in the "standard" setting (either dim=1,2 or 3), or the codimension 1 case (e.g. a 2-d manifold embedded in 3d space)? Or does one of your field variables only correspond to the codimension 1 case?

- Are you referring to some structures that store, for example, the kinematic variables for both fields as well as the tangents resulting from linearisation of the coupled problem?

- Are you solving your problem using a staggered scheme? I ask because you store neither any of the kinematic data nor tangents related to the thermal problem or coupling.

Have you examined the many tutorials that solve problem with multiple fields/components? They may not be multi-physics per se, but will give you some insights on how to set up your problem.

Regards,

J-P

[benhour....@gmail.com]

Dear All,

I want to solve a thermoelastic problem in cylindrical coordinate system. My domain is a circle in 2D or Sphere in 3D. I am following step 44 with considering the fact that I am solving for small strain. Because I have curved domain and I should Impose a boundary condition that change with temperature on the curved side of my domain, I understood that I should use MappingQ because of the fact that my whole domain includes curved edge. For this case, I have some problem with defining elasticity matrices such as thermal, initial and elastic strain along with stresses. I have attached my code for your consideration. It would be very kind of you if you take a glance at it and tell me what is wrong defining my matrices.

[Jean-Paul Pelteret]

Dear Benhour,

Thank you for providing a little more context, but you still have not articulated what the actual problem is. What exactly do you mean by "I have some problem with defining elasticity matrices"?

Please see this post, in particular the section titled "How to best formulate a question".

Regards,

J-P

[benhour....@gmail.com]

Dear Jean-Paul,

Let me put it in a such way. I want to solve my equation on a curved edge. According to dealii library, I should use MappingQ whenever I have nonlinear domain. Am I correct? In addition, Should I replace all dim in defining my tensors to spacedim? It should be noted that dim = spacedim - 1.

Best,

Benhour 

[Timo Heister]

> Let me put it in a such way. I want to solve my equation on a curved edge.
> According to dealii library, I should use MappingQ whenever I have nonlinear
> domain. Am I correct?

I don't think "nonlinear" is the word you are looking for. Do you mean
a curved boundary?

> In addition, Should I replace all dim in defining my
> tensors to spacedim? It should be noted that dim = spacedim - 1.

Using a mapping (see step-10, step-11, etc.) has nothing to do with
codim 1 or 2 problems (where dim is not equal to spacedim, see
step-38).

--
Timo Heister
http://www.math.clemson.edu/~heister/

[benhour....@gmail.com]

Yes You are right. I mean curved domain. My question is should i use codim instead of mapping because of the fact that there is only two tutorials in dealii which solved a differential equation on c curved domain according to the fact that my derived weak expressions are in cylindrical coordinate systems?

Thanks,

Benhour

[Timo Heister]

> Yes You are right. I mean curved domain. My question is should i use codim
> instead of mapping because of the fact that there is only two tutorials in
> dealii which solved a differential equation on c curved domain according to

No, codim problems solve systems on surfaces embedded in a higher
dimensional space. This has nothing to do with solving on a domain
that is not a rectangle. There are other examples that use a different
domain, see step-32 for example.

> the fact that my derived weak expressions are in cylindrical coordinate
> systems?

Are you sure you want to describe your problem in cylindrical
coordinates? While possible, this probably only makes sense if you
have additional symmetries in your system that you want to exploit,
but you will have many other issues with that (quadrature etc.).

If you just want to solve on a disk or sphere, you just need to change
the geometry and potentially use a higher order mapping for accuracy.

[benhour....@gmail.com]

Yes exactly. The only difference is in weak form of my equation which should be multiplied by 2*pi*r which I define r_value = fe_values_eta.quadrature_point(q_point)[0];

In fact, I should solve my problem on a sphere (circle). According to your noteworthy comments, the only thing that I have to do is change the geometry and use higher order mapping for accuracy. Another thing is that consider a half of a circle in r-z coordinate system and the problem is axisymmetric. what should be my boundary condition for rotation axis (z-coordinate)? For the r-coordinate, I can use no-normal flux boundary conditions or just limit the displacement normal to the surface, however, for rotation axis I am really confused. 

[Wolfgang Bangerth]

> In fact, I should solve my problem on a sphere (circle). According to your
> noteworthy comments, the only thing that I have to do is change the geometry
> and use higher order mapping for accuracy. Another thing is that consider a
> half of a circle in r-z coordinate system and the problem is axisymmetric.
> what should be my boundary condition for rotation axis (z-coordinate)?

That depends on the equation, but if you were, for example, solving the
Laplace equation in cylindrical coordinates, then there would be no special
boundary condition along the z-axis (at r=0). That's because you will find
that when you do the integration by parts, you get boundary terms of the form

\int_{\partial \Omega} n.nabla u(r,z) varphi(r,z) 2*pi*r dr dz

But at the symmetry axis, r=0, so the integrand disappears at the boundary.

In other words, for most equations, there is nothing special you need to do at
the symmetry axis.

Best
W.

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

[benhour....@gmail.com]

Dear All,

I still have some problem in solving my equations on a curved domain (sphere in 3D or circle in 2D). I got error with this content:

invalid initialization of reference of type ‘dealii::Triangulation<2, 3>&’ from expression of type ‘dealii::Triangulation<3, 3>’

Does it mean that I should use different codim for my domain and Finite element spaces because of the fact that I should solve my problem and impose my boundary condition on a sphere.

Thanks,

Benhour

[Wolfgang Bangerth]

On 10/19/2016 08:42 PM, benhour....@gmail.com wrote:
> I still have some problem in solving my equations on a curved domain (sphere
> in 3D or circle in 2D). I got error with this content:
> invalid initialization of reference of type ‘dealii::Triangulation<2, 3>&’
> from expression of type ‘dealii::Triangulation<3, 3>’

You pass a 3d triangulation to a function that expects a
2d-triangulation-embedded-into-3d-space as argument.

Which function are you trying to call?

[benhour....@gmail.com]

Dear All,

I want to solve my thermoelasticity problem in cylindrical coordinate systems. my domain includes two quarter of a circle in 2D and two quarter of two spheres in 3D. They are bulk. I attached my domain that was created by quarter_hyper_shell for 2D. But, I really do not kno how I can create it in 3D. It would be very kind of you if you help me to create my geometry according to the fact that the hyper_sphere only sxists in dim=spacedim + 1, so by considering this assumption, what changes should i impose on my equation? Should I change my DoFHandler or Trian gulation as well?

Bests,

Benhour  

[Valentin Zingan]

Hello Wolfgang,

Instead of 2*pi*r  dr  dz    in the boundary integral, did you mean 2*pi*r*dz ?

Otherwise, you will get the differential volume.

Thanks. 

среда, 19 октября 2016 г., 17:43:59 UTC-4 пользователь Wolfgang Bangerth написал:

[Wolfgang Bangerth]

> Instead of 2*pi*r  dr  dz    in the boundary integral, did you mean 2*pi*r*dz ?
> Otherwise, you will get the differential volume.

This is a really old thread :-) Yes, you are correct. The *domain* integral
has the form

\int_\Omega f(r,z) 2*pi r dr dz

and the boundary integral will have the form

\int_{\partial\Omega} f(R,z) 2*pi R dz

Best
W.

> среда, 19 октября 2016 г., 17:43:59 UTC-4 пользователь Wolfgang Bangerth написал:
>
> > In fact, I should solve my problem on a sphere (circle). According to your
> > noteworthy comments, the only thing that I have to do is change the
> geometry
> > and use higher order mapping for accuracy. Another thing is that
> consider a
> > half of a circle in r-z coordinate system and the problem is axisymmetric.
> > what should be my boundary condition for rotation axis (z-coordinate)?
>
> That depends on the equation, but if you were, for example, solving the
> Laplace equation in cylindrical coordinates, then there would be no special
> boundary condition along the z-axis (at r=0). That's because you will find
> that when you do the integration by parts, you get boundary terms of the form
>
>    \int_{\partial \Omega}   n.nabla u(r,z)  varphi(r,z)  2*pi*r  dr  dz