The answer doesn't converge by refining the mesh

General Discussion
1

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"
> mesh Th=gmshload("Chiral.msh");
> 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?
Thank you in advance for your answers.

2

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.

3

Thank you for your answer.
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))

and the sigma11 had converged.
thanks again.
FreeFem++.rar (208.5 KB)

4

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

5

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"
//mesh Th=gmshload("Chiral.msh");
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)

6

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;