# The answer doesn't converge by refining the mesh

Hi, I am trying to solve a BVP to obtain the micropolar strain localization tensor in a 2D chiral geometry with some periodic boundary conditions. I have 2 unknowns h111 and h211.
The boundary conditions are:
h111(-0.5,y)= h111(0.5,y)
h211(-0.5,y)= h211(0.5,y)
h111(x,-0.5)= h111(x,0.5)
h211(x,-0.5)= h211(x,0.5)
Below is the code I’ve been using:

``````> load "gmsh"
> plot(Th,wait=1);
> fespace Vh(Th,P1,periodic=[[18,x],[20,x],[19,y],[21,y]]);
> Vh h211,w2,h111,w1,b=24500/13,c=10500/13,d=14000/13;
> solve diffuse([h111,h211],[w1,w2])=
> int2d(Th)(dx(w2)*(d*dx(h211))
> + dx(w1)*(b*dx(h111) + c*dy(h211))
> + dy(w1)*(d*dx(h211))
> + dy(w2)*(c*dx(h111) + b*dy(h211)))
> + int2d(Th) (b*dx(w1))+int2d(Th) (c*dy(w2));
> real sigma11=int2d(Th)(((24500*dx(h111))/13 + 24500/13)*(dx(h111) + 1)
> + (10500*dy(h211)*(dx(h111) + 1))/13 + (28000*dx(h211)^2)/13 + (24500*dy(h211)^2)/13 + dy(h211)*((10500*dx(h111))/13 + 10500/13));
>
> cout << "sigma11 " << sigma11 << endl;
``````

However, sigma11, calculated using h111 and h211, does not converge by refining the mesh. Sometimes it decreases and sometimes it increases. What should I do?

Can you post the code as an attachment as my browser makes copy paste
screw ups. I don’t know if I’ll look at it but curious. Do you have a link to areference
or the eqn you are trying to solve? Probably it does not matter a lot on the mesh
for diagnosis although if you have specific features of the mesh it may help to post it.
Thanks.

The code and the geometry are attached to the post.
In my equation, dx(h211) equals dy(h111) because the micropolar strain localization tensor is symmetric. Before I considered that this tensor is symmetric, my equation was:

`````` int2d(Th)(dx(w2)*(a*dx(h211) + a*dy(h111))
+ dx(w1)*(b*dx(h111) + c*dy(h211))
+ dy(w1)*(a*dx(h211) + a*dy(h111))
+ dy(w2)*(c*dx(h111) + b*dy(h211)))
+ int2d(Th) (b*dx(w1))+int2d(Th) (c*dy(w2))
``````

thanks again.
FreeFem++.rar (208.5 KB)

The geometry.geo is attached to this post.
Chiralgeo.rar (729 Bytes)

I don’t have time to sort it out as there are a lot of constants in the signal expression
and you didn’t post the equation. I get a constant signal over a range of mesh values
using a simple square mesh,
See if the attached code can help you sort it out. Also did you intend to have a
cross product in there- could be a sign mistake.

``````//load "gmsh"
int nn=5;
mesh Th = square(nn,nn,[100*x,100*y]);
//fespace Vh(Th, P1); //scalar FE
fespace Vh(Th,P2,periodic=[[2, y],[4, y],[1, x],[3, x]]); //bi-periodic FE
//fespace Vh(Th,P2,periodic=[[18,x],[20,x],[19,y],[21,y]]);
Vh h211,w2,h111,w1,a=7000/13,b=24500/13,c=10500/13,d=14000/13;
solve diffuse([h111,h211],[w1,w2])=
int2d(Th)(dx(w2)*(d*dx(h211))
+ dx(w1)*(b*dx(h111) + c*dy(h211))
+ dy(w1)*(d*dx(h211))
+ dy(w2)*(c*dx(h111) + b*dy(h211)))
+ int2d(Th) (b*dx(w1))+int2d(Th) (c*dy(w2));
real sigma11=int2d(Th)(((24500*dx(h111))/13 + 24500/13)*(dx(h111) + 1)
+ (10500*dy(h211)*(dx(h111) + 1))/13 + (28000*dx(h211)^2)/13 + (24500*dy(h211)^2)/13 + dy(h211)*((10500*dx(h111))/13 + 10500/13));

plot(h111,wait=1,fill=1);
plot(h211,wait=1,fill=1);
cout << "sigma11 " << sigma11 << endl;

``````

test4.edp (853 Bytes)

Also I would get in the habit of writing expressions with decimal points
as the consts may be evaluated as int/int but probably that is just a small
issue. You can check with this though,

``````cout<<"c="<<c<<endl;

``````