I want to implement elasticity in domain [0,1]*[0,1] with known analytical solution, and I want to use different values of lambda and nu, in different halves of the domain and I tried it as follows:
real nu1 = 10;
real lambda1 = 2*10^4; Lame constant
real nu2 = 20;
real lambda2 = 10^4;
func nu = (y>=0.5 )*(nu1 + 0*x + 0*y) + (y<0.5 )*(nu2+ 0*x + 0*y);
func lambda = (y>=0.5 )*(lambda1 + 0*x + 0*y) + (y<0.5 )*(lambda2+ 0*x + 0*y);
On doing this, I noticed that at y=0.5, the L2 norm of computed solution u and exact solution uex are very different and I do not get the convergence of my solution. When I use the same parameters nu and lambda, I do not have any problem with the elasticity solver (I get desired convergence rates).
Please let me know if I needed to change anything in defining nu and lambda differently in different parts of the domain.
In the documentation section 2.6 (thermal conduction problem), different diffusion (k) has been used in different parts of the domain. So, I thought it should work the same way!
1)Defining nu and lambda as functions of x and y.
2) Use it in defining right-side functions and in variational formulation.
I plotted my solution in the case of elasticity in the whole domain [0,1]*[0,1] and it is way off as compared to exact solution. While using same lambda and nu over entire domain, my solution matches with the exact solution.
//RHS for Elasticity
func f1s = -nus*(2*u1struexx+u1strueyy+u2struexy)-lambdas*(u1struexx+u2strueyx);
func f2s = -nus*(2*u2strueyy+u2struexx+u1strueyx)-lambdas*(u1struexy+u2strueyy);
fespace Uh(Gh,P2); // element space for displacement
Uh u1s,u2s, v1s, v2s;
problem structure([u1s,u2s], [v1s,v2s]) =
int2d(Gh) (
2*nus*(
dx(u1s)*dx(v1s)
+ 0.5*(dx(u2s) + dy(u1s))*(dx(v2s) + dy(v1s))
+ dy(u2s)*dy(v2s)
) // 2nu(D(u),D(v))
+ lambdas*(dx(u1s) + dy(u2s))*(dx(v1s) + dy(v2s)) // lambda(div(u),div(v))
)
- int2d(Gh)(
f1s*v1s + f2s*v2s
)
+ on(2,3,4,5,u1s=u1strue,u2s=u2strue); // BC
Your “manufactured solution” does not solve the problem you want to solve, because of the jumps in \nu and \lambda. Indeed one has -\mathop{\rm div}(2\nu Du)-\nabla(\lambda\mathop{\rm div} u) =-2\nu\mathop{\rm div} Du-\lambda\nabla\mathop{\rm div} u - 2Du\nabla\nu-(\mathop{\rm div }u)\nabla\lambda and \nabla\nu=\delta(y-1/2)[\nu](0,1), \nabla\lambda=\delta(y-1/2)[\lambda](0,1),
where [\nu] denotes the jump of \nu through the line y=1/2, \delta(y-1/2) is the Dirac function on the line, and (0,1) is the unit vector oriented upwards.
Thus to be correct you must add a right-hand side \delta(y-1/2)g.
This is what I have done with the functions g1s, g2s, and defining the line as border G5 elasticdisc.edp (3.0 KB)
Thank you for your time and response @fb77. It really means a lot to me.
I was wondering why the domain-dependent parameters in the documentation are NOT handled similarly but they work (thermal problems). I didn’t see any right-side adjustment in those problems. And also, in this simple problem too, the errors look good now but interestingly when I plot the solution to compare with exact function, they are very OFF, even the scale, (see the attached plots for x-components):
@fb77, the way line y=1/2 was created in the mesh overlapped the mesh twice, so the solution looked weird. I modified it as mentioned above by creating Gh1 and Lh and glued them to create Gh. The solution is now very similar to the exact solution. THANK YOU SO MUCH FOR HELPING ME WITH THIS.
the domain-dependent parameters in the documentation are NOT handled similarly but they work (thermal problems).
The variational formulation is done here with nothing special. Only the data (right-hand side) is special. This special right-hand side is chosen in order to have the exact solution. In the thermal problems of the documentation, there is no exact solution. You can take an arbitrary right-hand side (with or without line integral) and get the corresponding solution numerically, this what is done in the examples of the documentation.
Another way of building an exact solution which is more “natural”, without line integral in the right-hand side, would be for the elliptic problem -\mathop{\rm div}(k\nabla u) =f, with k discontinuous, to take u_{ex} continuous, but with \nabla u_{ex} discontinuous in such a way that k\nabla u_{ex} is continuous. Then the corresponding f would have no line Dirac part.
Indeed when -\mathop{\rm div}(k\nabla u) =f with k discontinuous, u and f cannot be both smooth; at least one of them must have some kind of unsmoothness.
