Re: [Getfem-users] Displacements are mesh-size dependent in a simple 2D linear elastic problem?

[Yves Renard]

Dear Tom,

My opinion is that the definition of your pressure force density is not correct 
because you multiply the force density by the length of the tangent and because 
you integrate on the whole domain and not on the boundary. You should pass to  
getfem::add_source_term_brick the contact pressure (you can use also directly 
getfem::add_normal_source_term_brick in your case) and give the as an 
additional argument the boundary number on which you want to prescribe your 
Neumann condition (as you did for the Dirichlet condition).

Yves.

----- Original Message -----
From: "Tom Haeck" <address@hidden>
To: "Yves Renard" <address@hidden>
Sent: Friday, December 19, 2014 10:16:18 AM
Subject: Re: [Getfem-users] Displacements are mesh-size dependent in a simple 
2D linear elastic problem?

Hey Yves,

Thanks for your response!  As for now, I was calculating the length based on 
the actual physical point locations.  Although this does the job, I will try 
the more general way you suggested.

Concerning the displacements, I meant that given the forces on the faces and 
given the appropriate boundary conditions, the linear elastic FEM model is 
solved for the displacements.  I include a small code snippet to clarify:

    getfem::mesh_im mim(mesh);
    getfem::mesh_fem mf_u(mesh);
    getfem::mesh_fem mf_rhs(mesh);



    scalar_type lambda, mu;

    mf_u.set_qdim(dim_type(N));



    mim.set_integration_method(getfem::int_method_descriptor("IM_TRIANGLE(6)"));



    mf_u.set_finite_element(getfem::fem_descriptor("FEM_PK(2,1)"));
    mf_rhs.set_finite_element(getfem::fem_descriptor("FEM_PK(2,1)"));



    // solve
    plain_vector U;
    getfem::model model;



    // Main unknown of the problem.
    model.add_fem_variable("u", mf_u);



    // Linearized elasticity brick.
    model.add_initialized_scalar_data("lambda", lambda);
    model.add_initialized_scalar_data("mu", mu);
    getfem::add_isotropic_linearized_elasticity_brick(model, mim, "u", 
"lambda", "mu");



    // Volumic source term.
    std::vector<scalar_type> F(mf_rhs.nb_dof()*N);
    ComputeForces(F,mesh,mf_rhs,pressure);
    model.add_initialized_fem_data("VolumicData", mf_rhs, F);
    getfem::add_source_term_brick(model, mim, "u", "VolumicData");



    // Dirichlet condition.
    gmm::resize(F, mf_rhs.nb_dof()*N);
    for(int i=0; i<F.size(); i++) F[i]=0.0;
    model.add_initialized_fem_data("DirichletData", mf_rhs, F);
    getfem::add_Dirichlet_condition_with_multipliers(model, mim, "u", mf_u, 
OUTER_BOUNDARY, "DirichletData");



    // Solve
    gmm::iteration iter(residual, 1, size_type(-1));
    iter.init();
    getfem::standard_solve(model, iter);
    gmm::resize(U, mf_u.nb_dof());
    gmm::copy(model.real_variable("u"), U);
    if (!iter.converged()) GMM_ASSERT1(false, "Solve has failed. No 
convergence.”);


Corresponding to a uniform pressure, the forces are calculated in the separate 
function “ComputeForces”:

void ComputeForces(plain_vector& forces, getfem::mesh& mesh, getfem::mesh_fem& 
mf_rhs, scalar_type pressure)
{
    ...
    for (getfem::mr_visitor i(mesh.region(INNER_BOUNDARY)); !i.finished(); ++i) 
{
        bgeot::base_node dir = mesh.normal_of_face_of_convex(i.cv(), i.f());
        scalar_type dirMag = gmm::vect_norm2(dir);
        dir /= -dirMag;



        bgeot::base_node tangent(N);
        tangent[0] = (mesh.points_of_face_of_convex(i.cv(), i.f())[0][0] -  
mesh.points_of_face_of_convex(i.cv(), i.f())[1][0]);
        tangent[1] = (mesh.points_of_face_of_convex(i.cv(), i.f())[0][1] -  
mesh.points_of_face_of_convex(i.cv(), i.f())[1][1]);
        scalar_type tanMag = gmm::vect_norm2(tangent);
        dir *= (tanMag*pressure);
        dir /= (scalar_type) N;



        for (int j = 0; j < N; j++)
            for (int k = 0; k < N; k++)
                forces[mf_rhs.ind_basic_dof_of_face_of_element(i.cv(), 
i.f())[j]*N + k] += dir[k];
    }
    return;
}

For this code snippet, the solution U is heavily dependent on the sizes of the 
elements in the mesh, which makes me think that I made some mistake…

Moreover, I now try to introduce a generalized Dirichlet boundary condition, 
H*u=r:

plain_vector H(mf_rhs.nb_dof()*4); // fill with values
plain_vector r(mf_rhs.nb_dof()*2); // fill with values
model.add_initialized_fem_data(“Hdata",mf_rhs, H);
model.add_initialized_fem_data(“rdata",mf_rhs, r);

getfem::add_generalized_Dirichlet_condition_with_multipliers(model, mim, "u", 
mf_u, OUTER_BOUNDARY , “Hdata", “rdata");

This gives an error:
libc++abi.dylib: terminating with uncaught exception of type gmm::gmm_error: 
Error in getfem_assembling_tensors.cc, line 811 :
tensor error: wrong number of indices for the 3th argument of the reduction 
comp(vBase(#1).vBase(#3).Base(#2)(:,i,:,j,k).F(i,j,k)) (expected 1 indexes, got 
5

I’d really appreciate it if you could take a look at my code.

Greetings,

Tom


On 16 Dec 2014, at 17:42, Yves Renard <address@hidden<mailto:address@hidden>> 
wrote:


Deat Tom,

Yes,  normal_of_face_of_convex gives a normal vector which lenght is the 
determinant of the transformation from the reference face on the real face. 
However ... the area of the face of the reference element depends if it if a 
face parallel to an axis (lenght 1) or the third face wich is of length 
sqrt(2). So the length of the normal vector is not directly the area (length) 
of the face.

If I understand well, you want to compute the equivalent forces on each face 
corresponding to your uniform pressure.

One simple way to obtain the lenght/area of a face (even for curved ones) is to 
sum the weights of an integration method on the corresponding face multiplied 
by the determinant of the transformation (i.e. to integrate 1 on the face). If 
you do not want to go to the internal of getfem, you can simply ask to compute 
the mass matrix between a P0 finite element method (eventually reduced to the 
boundary you look at) and a finite element on which you define the pressure. Or 
if you use the high level generic assembly, you can simply ask the assembly of 
the term "-Normal*Pressure*Test_p0" on the corresponding boundary where 
Pressure is your contact pressure and p0 is a field defined on a P0 finite 
element method.

Concerning the displacements, what do you mean by "I calculate the resulting 
displacement ?". You calculate some displacements explicitely from the computed 
force on each face ? I do not see very well how you can do that. Could you 
explain more this point ?

Yves.



Le 15/12/2014 16:32, Tom Haeck a écrit :
Hi all,

I am experimenting with a very simple 2D linear elastic problem.  The domain is 
(more or less) rectangular with a circular hole in the middle.  A constant 
pressure is acting radially outward on the faces of the hole in the middle (see 
attached image).

In order to calculate the force that is acting on each face of the hole in the 
middle, I multiply the constant pressure by the length of the face.  Next, I 
calculate the resulting displacements.  It struck me that the magnitude of the 
resulting displacement is largely dependent on the size of the elements of the 
mesh.  I would expect that larger forces acting on a few larger mesh elements 
would give more or less the same displacements as smaller forces acting on a 
lot of small mesh elements?

Moreover, when I use the magnitude of the vectors 
mesh.normal_of_face_of_convex(i.cv(),i.f()) as magnitude of the force vectors 
acting on the faces of the hole, the resulting displacements are independent of 
the mesh size.  The magnitudes of the vectors mesh.normal_of_face_of_convex() 
are larger for smaller elements and smaller for bigger elements.  However, 
according to the documentation,  mesh.normal_of_face_of_convex() returns 
vectors with magnitudes that are proportional to the area of the face?

To summarize:

  *   how come that my displacements are dependent on the element mesh size?
  *   how come that the displacement are not dependent on the element mesh size 
anymore when I use normal_of_face_of_convex?
  *   how come that the magnitudes of the vectors of normal_of_face_of_convex 
are not proportional to the area of the face, as is documented?

Thx,

Tom


<Mail Attachment.png>




_______________________________________________
Getfem-users mailing list
address@hidden<mailto:address@hidden>
https://mail.gna.org/listinfo/getfem-users




--

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

---------