# Mesh problem / Auto contact

**URL:** <https://community.freefem.org/t/mesh-problem-auto-contact/3146>\
**Category:** General Discussion\
**Created:** [April 24, 2024, 3:46pm UTC](https://community.freefem.org/t/mesh-problem-auto-contact/3146 "2024-04-24T15:46:12Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![Joujou](https://avatars.discourse-cdn.com/v4/letter/j/5e9695/32.png) [@Joujou](https://community.freefem.org/u/Joujou)\
**Post date:** [April 24, 2024, 3:46pm UTC](https://community.freefem.org/t/mesh-problem-auto-contact/3146/1 "2024-04-24T15:46:13Z")

</div>

Hello everyone,

I am working on a simulation regarding the compression of a pillar ( Rectangular form).

At high level of deformation, two elements touch each other and my calculation get stopped(see the figure attached).

Please any suggestion is very welcome.

![pastedImage1](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/9/94362090c4be5d46fce715a12f0ba34d14797442.png)

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [April 24, 2024, 4:45pm UTC](https://community.freefem.org/t/mesh-problem-auto-contact/3146/2 "2024-04-24T16:45:42Z")

</div>

Hello Jalal,  
It is not clear what your picture means. What do the colors represent? I understand that the mesh is deformed and that your problem is its boundary.  
We had a similar problem in  
J. Comput. Phys., 333:387–408, 2017. DOI: 10.1016/j.jcp.2016.12.036  
[https://perso.math.u-pem.fr/bouchut.francois/publications/drucker-prager2D.pdf](https://perso.math.u-pem.fr/bouchut.francois/publications/drucker-prager2D.pdf)  
We managed to solve it by putting some numerical surface tension, see in particular figure 4 p. 11, to prevent the folding of the free surface over itself.  
This surface tension is the last term in the right-hand side of the variational formulation 2.17a

---

<div class="post-metadata">

**Author:** ![Joujou](https://avatars.discourse-cdn.com/v4/letter/j/5e9695/32.png) [@Joujou](https://community.freefem.org/u/Joujou)\
**Post date:** [April 25, 2024, 6:31am UTC](https://community.freefem.org/t/mesh-problem-auto-contact/3146/3 "2024-04-25T06:31:39Z")

</div>

Thank you François for your help.

I found the paper very interesting and aligned perfectly with what I was searching for.  
I am trying to implement the numerical tension.

The simulation is a compression test of a multilayer pillar, made of alternation of two compositions.

The colors in the picture (Blue and Orange) correspond to the yield stength values of the two materials (Compositions). I use this representation to identify the layers during the simulation animation’s.

 ![pastedImage](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/9/9d8f294c099b0d91a93881b8035e7dd39420178f.png)

It is a coarse mesh indeed but I have the same problem when the mesh is fine.

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [April 25, 2024, 12:34pm UTC](https://community.freefem.org/t/mesh-problem-auto-contact/3146/4 "2024-04-25T12:34:36Z")

</div>

You make a nice computation!  
In our code I remember that obtaining the local curvature of the boundary was a bit messy.  
Here is the piece of code (sorry comments are in french). It is not optimal:  
the construction of R is repeated on 3 pieces of the boundary

```auto
  real[int] rx(Th.nbe),ry(Th.nbe),curvx(Th.nbe), curvy(Th.nbe);
  int[int] ibe(Th.nbe);
  real curvx0,curvy0;
  int nbv=0; // compte le nombre de noeud au bord
  // les noeuds du bord sont ranges dans l'ordre: label=1,2,3,4 
  // surface libre parcourue de droite (2) à gauche (3)
  for(int n=0;n<Th.nbe;n++){
    if(Th.be(n).label==2){
      //-------------------- construction R(x,y) ----------------------
      rx(nbv)=Th.be(n)[0].x-Th.be(n)[1].x;
      ry(nbv)=Th.be(n)[0].y-Th.be(n)[1].y;
      curvx(nbv)=Th.be(n)[1].x;
      curvy(nbv)=Th.be(n)[1].y;
      ibe(nbv)=n;
      if(nbv==0){
        curvx0=Th.be(n)[0].x;
	curvy0=Th.be(n)[0].y;
      }
      nbv++;
    }
  }
  int ncrit=nbv-1;
  for(int n=0;n<Th.nbe;n++){
    if(Th.be(n).label==31){
      //-------------------- construction R(x,y) ----------------------
      rx(nbv)=Th.be(n)[0].x-Th.be(n)[1].x;
      ry(nbv)=Th.be(n)[0].y-Th.be(n)[1].y;
      curvx(nbv)=Th.be(n)[1].x;
      curvy(nbv)=Th.be(n)[1].y;
      ibe(nbv)=n;
      if(nbv==0){
        curvx0=Th.be(n)[0].x;
	curvy0=Th.be(n)[0].y;
      }
      nbv++;
    }
  }
  for(int n=0;n<Th.nbe;n++){
    if(Th.be(n).label==32){
      //-------------------- construction R(x,y) ----------------------
      rx(nbv)=Th.be(n)[0].x-Th.be(n)[1].x;
      ry(nbv)=Th.be(n)[0].y-Th.be(n)[1].y;
      curvx(nbv)=Th.be(n)[1].x;
      curvy(nbv)=Th.be(n)[1].y;
      ibe(nbv)=n;
      if(nbv==0){
        curvx0=Th.be(n)[0].x;
	curvy0=Th.be(n)[0].y;
      }
      nbv++;
    }
  }
  rx.resize(nbv); ry.resize(nbv);
  curvx.resize(nbv); curvy.resize(nbv);
  ibe.resize(nbv);
  real[int] arctan(nbv+1),curv(nbv+1),num(nbv+1),denom(nbv+1),lg(nbv+1);
  real normr,normm;
  curv(0)=0;
  for(int n=0;n<nbv-1;n++){
    num(n+1)=-ry(n+1)*rx(n)+rx(n+1)*ry(n); // rt(n).r(n-1) avec rt(n)=(-ry,rx)
    normr=sqrt(rx(n+1)^2+ry(n+1)^2);
    normm=sqrt(rx(n)^2+ry(n)^2);
    denom(n+1)=normm*normr+rx(n+1)*rx(n)+ry(n+1)*ry(n);
    arctan(n+1)=2*atan(num(n+1)/denom(n+1));
    lg(n+1)=0.5*(normr+normm);
    curv(n+1)=arctan(n+1)/lg(n+1);
    //cout << "curv("<<n<<")="<<curv(n)<<endl;
  }
  curv(nbv)=curv(nbv-1);
  //-------------------- COURBURE ---------------------
  real h=0.025, cst=0., deltan=16.;
  fcurv[]=0.; fcurv=0.;
  fcurv[][int(Th.be(ibe(0))[0])]=curv(0);
  for(int n=0;n<nbv;n++)
    fcurv[][int(Th.be(ibe(n))[1])]=h*curv(n+1)*max(0.,1.-((n-ncrit)/deltan)^2)^2; 

  {
    ofstream fichier("curvature.dat");
    for(int n=0;n<nbv;n++){
      fichier << curv(n) << " " << ibe(n) << " " << curvx(n) << " " << curvy(n) << " " << num(n) << " " << denom(n)<< " " << arctan(n)<< " " << lg(n)<<"\n"; 
    }
  }

```

full code  
[drucker-prager\_ale.edp](https://community.freefem.org/uploads/short-url/5FT0VfJCLohJj82fKXaBqfXTeWK.edp) (9.9 KB)