Re: [Getfem-users] problem with structural calculations in getfem

[Yves Renard]

On Tuesday 06 November 2007 17:50, Andriy Andreykiv wrote:
> Dear Yves,
>
>     First of all thank you very much for debugging that two NonLin term
> problem. If I also may consult you regarding the following:
>
>     While looking for bugs in my electrostatic-structural implementation I
> confrontated a problem with quite simple structural calculation in getfem.
> Although this calculation gives a result that looks nice, comparing the
> numbers whith commercial FEM package MSC. Marc gave very different
> solution. The problem is very simple (please see the picture in the
> attachment):
>
> A rectange with a size 0.8 by 0.1 is clamped at the sides.
> A constant distributed load P is applied on the bottom surface of this
> rectangle. In getfem, this load (pressure) is prescribed by source_term,
> with  vector "V(#1)+=comp(Normal().vBase(#1))(i,:,i);", assembled on the
> bottom boundary. Then this vector is scaled with the actual value of the
> pressure.  In both Marc and Getfem the corresponding nodal force is
> calculated on the underformed geometry of the body (hence, no follower
> force). The force causes the rectangle to subside downwards and the
> deflection of the middle part of the rectangle is measured.
>
> Elastic properties (used by MSC. Marc) are
>                                Young moduls: 20
>                                 Poisson ratio: 0.3
>                                 Plain strain assumed for 2D problem
>
> Corresponding Lame constants (used in Getfem program):
>                                 Lambda:    11.5384
>                                        Mu:    7.6923
>
> In both programs I used linear quad elements, splitting the rectangle 32X4.
> Both programs use 4 Gauss points for element integration. A similar problem
> was also solved in 3D for 0.8X0.1X0.1 paralellepiped structure, descretized
> by 32X4X4 linear hex elements. Bellow you can see small and finite strain
> results from the two programs:
>
> 2D, finite strain, distributed load P=-2.0:
>         MSC. Marc deflection = 0.229; Getfem deflection 0.192651
>
> 2D, small strain, distributed load P=-0.4
>         MSC.Marc deflection = 0.2667; Getfem deflection 0.1699
>
> 3D, finite strain, distributed load P=-2.0
>         MSC. Marc deflection = 0.2375; Getfem deflection 0.0191458 !!!
>
> 3D, small strain, distributed load P=-0.4
>         MSC. Marc deflection = 0.287; Getfem delflection 0.00400608 !!!
>
> Even if Marc results are not correct, the agreement between 2D and 3D
> solution is much better then in Getfem. Additionally, when doing either p-
> or h- refinement in Marc, the result changes only slightly, by increasing
> the deflection (2-3 %). When I refine the mesh in Getfem, the deflection
> considerably reduces which is really strange.
>
> The source of the program with the corresponding four parameter files
> (2D/3D and finite/small strain cases) can be found in the attachment.
>
> Thank you in advance,
>                                     Andriy
>
>
>
>
> __________________________________________
> Andriy Andreykiv (PhD, MSc)
> Delft University of Technology
> Faculty of Mechanical Engineering
> Material Science and Marine Technology
> Group: Fundamentals of Microsystems
> Mekelweg 2
> 2628 CD  Delft
> The Netherlands
>
> E-mail  :  address@hidden
> Tel.    :   (31) 15 2786818
> Fax.    :   (31) 15 2789475
> www     :   http://www-tm.wbmt.tudelft.nl/~andrico
> private:  (31) 6 47376804


Dear Andriy,

I see at least one error in your program. When you run the following lines

plain_vector Unit_Pressure(mf_mechanic.nb_dof());
getfem::generic_assembly assem_boundary;
assem_boundary.set("V(#1)+=comp(Normal().vBase(#1))(i,:,i);");
assem_boundary.push_mi(mim_mechanic);
assem_boundary.push_mf(mf_mechanic);
assem_boundary.push_vec(Unit_Pressure);
assem_boundary.assembly(m_mechanic.region(SOUTH));

plain_vector Current_Pressure(Unit_Pressure);
gmm::scale(Current_Pressure,p_current);
pstruct_brick->source_term().set(Current_Pressure);

The vector you give to pstruct_brick->source_term().set(.) is already 
assembled. But, it will be assembled a second time ! You dont have to give an 
assembled vector to this function. You have to give a vector of degree of 
freedom corresponding to your mf_data. For instance a more correct code will 
be

plain_vector Current_Pressure(data_size*m_mechanic.dim());
gmm::fill(gmm::sub_vector(Current_Pressure,
                  gmm::sub_slice(1, data_size, 2)),
                  p_current);
pstruct_brick->source_term().set(Current_Pressure);

(the fact that the vector is defined outside the considered zone is not 
important).

Yves.

-- 

  Yves Renard (address@hidden)       tel : (33) 04.72.43.87.08
  Pole de Mathematiques,                       fax : (33) 04.72.43.85.29
  INSA de Lyon, Universite de Lyon
  20, rue Albert Einstein
  69621 Villeurbanne Cedex, FRANCE
  http://math.univ-lyon1.fr/~renard

---------