Linear elasticity 2D mesh: Error

Question

Hi,

The following is the code for solving elasticity (linear) problem on a simple square mesh. I took this from an old launchpad question.

https://github.com/ChaitanyaGoyal/2D_Linear_elasticity_Square_Plate/blob/master/Fenics_Drichilet_Load%20on%20top%20edge.py

It shows the following error when run:

page@Jarvis:~$ python 2D_Elasticity_node.py
Traceback (most recent call last):
File "2D_Elasticity_node.py", line 38, in
obstacle.mark(domains, 1)
File "/usr/lib/python2.7/dist-packages/dolfin/cpp/mesh.py", line 4623, in mark
self._mark(*args)
TypeError: SWIG director type mismatch expected 'bool' as the output argument of 'inside'

I can't see anything wrong with that statement as I have confirmed the domain coding from the following fenics documentation:

http://fenicsproject.org/documentation/dolfin/1.2.0/python/demo/pde/subdomains-poisson/python/documentation.html

Pls. help. Thanks

Answer 1

I think you forgot to actually put near(...) in the output of the inside() function for obstacle :D.

But then I get a different problem if I try to start your code because there is another error. In line 41, you set the "default" domain value to 0, but then in line 81 you integrate over domain 1. After replacing 1 by 0 everything works fine for me.

regards

Comments

Hi multigrid202, Thanks a lot for your input. The 'near' problem was correct.
After setting the default domain to 0, I set obstacle.domain to 1 and thats why I used dx(1)

This was suggested by Jan Blechta at the following link:

https://answers.launchpad.net/dolfin/+question/227859

However, Could you pls. explain when this guy says I want to apply a neumann BC on node (1,1) and then Jan Blechta asks him go for the dx(1) approach, what does this imply physically, what are they trying to do to node (1,1)?

Thanks a ton for taking interest and helping me out with these trivial conceptual doubts. Means a lot.

---

Hello there,

I don't know why they used dx(1) in that old question, it really doesn't make much sense. I suspect that it is a typo. The bilinear from $a$ is always defined on your domain and never on any obstacles. I also don't really get why the upper right corner is called obstacle ^^

Anyway, physically they are trying to apply a force to that corner. Imagine attaching a piece of string there and then pulling...

Two more remarks on that:

From a mathematical point of view, a pointwise boundary condition is meaningless. Physically, imagine pushing an infinitely thin needle into your desk. It will punch a hole at any non-zero force. This often is not a problem in a finite element code if the mesh is coarse enough, but personally, I consider it bad practice. If you have a fine mesh at the point where the force is applied, then you will get weird results.

The dx(1) has nothing to do with the boundary condition. dx(1) (or rather dx(0)) only specifies the domain where your problem is defined. But you only really need that if there is a reason to split your domain into two subdomains, e.g. if you have a plate that has a different youngs's modulus in one area. But in the current example, there is no real reason to use it. Just delete the definition of dx=... and only put dx where it says a=... Now, I explained in a different post that if you don't specify any boundary condition, then you will automatically get homogenous neumann conditions. If you want to apply a nonzero tractions force $g$ at a part of the boundary $\Gamma_N$, then you need to add $\int_{\Gamma_N}gv\;ds$ to the right hand side of your weak form ($v$ is a test function). Note that if $g=0$, this corresponds to "doing nothing" as explained above.

On another note: I strongly suggest that you take some type of course on the finite element method. Youtube is a good source. It has lectures e.g. from MIT. I am obviously happy to help, but I think that a structured introduction to the finite element method would be very beneficial for you. Especially, because you already know fenics and you can try out all the concepts during the lectures :)

best regards

---

Thank you very much for explaining that. I' ll post this link as a comment on the old launchpad question for anyone who might want to further pursue that thread in future.

I' ll definitely follow your suggestion the online learning for better understanding. Thanks!!