### Functionspace in 3D problem

377
views
0
7 months ago by
Hello,
I work on a 3D elasticity project. I want to create a 3D functionnal space with V, but when the gradient of a function in this space seems to be a (3,2) dimension matrix, so I deduced the the function is a 2D vector. How the define a 3D functional space (it doesn't seems to work in my code ) ?
my code :
[...]
V = VectorFunctionSpace(mesh,"P",1,dim=3)
u_ = TestFunction(V)
v = TrialFunction(V)
u = Function(V, name="Deplacement")

def eps(w):
def sigma(w):

return lmbda*tr(eps(w))*Identity(2) + 2*mu*eps(w)

Wdef = inner(sigma(v), eps(u_))*dx
Wext = dot (ps, u_)*ds(6)​
and the error message:
Cannot take symmetric part of rectangular matrix with dimensions (3, 2).
Traceback (most recent call last):
File "script_tube.py", line 58, in <module>
Wdef = inner(sigma(v), eps(u_))*dx
File "script_tube.py", line 55, in sigma
return lmbda*tr(eps(w))*Identity(2) + 2*mu*eps(w)
File "script_tube.py", line 52, in eps
File "/usr/lib/python2.7/dist-packages/ufl/operators.py", line 294, in sym
return Sym(A)
File "/usr/lib/python2.7/dist-packages/ufl/tensoralgebra.py", line 466, in __new__
error("Cannot take symmetric part of rectangular matrix with dimensions %s." % (sh,))
File "/usr/lib/python2.7/dist-packages/ufl/log.py", line 172, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Cannot take symmetric part of rectangular matrix with dimensions (3, 2).
Aborted (core dumped)​
Best Regards and thank you.
Community: FEniCS Project
2
Is your domain 2D or 3D?
If you're working in 2D, your displacement should also have dim=2. Or just omit the dim kwarg in the VectorFunctionSpace creation and let it derive the dim from your mesh argument.

The gradient of a 2D vector on a 2D domain is a 2x2 matrix, the gradient of a 3D vector on a 3D domain is a 3x3 matrix. The gradient of a 3D vector on a 2D domain is a 3x2 matrix.
written 7 months ago by klunkean
ok I understand, I modified my code an the error is gone. But an other dimension problem appear with the function dot(.,.)
I defined the vector ps like
ps = Fs*Constant((0,-1,0))
​
and in the command (u_ is defined in my previous code )
Wext = dot (ps, u_)*ds(6)
​
the following error message is raised :
 File "script_tube.py", line 60, in <module>
Wext = dot (ps, u_)*ds(6)
File "/usr/lib/python2.7/dist-packages/ufl/operators.py", line 163, in dot
return Dot(a, b)
File "/usr/lib/python2.7/dist-packages/ufl/tensoralgebra.py", line 213, in __new__
error("Dimension mismatch in dot product.")
File "/usr/lib/python2.7/dist-packages/ufl/log.py", line 172, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Dimension mismatch in dot product.
Aborted (core dumped)​
written 7 months ago by Liam Gonzalez
You should be able to check your dimensions by using the ufl_shape() method on objects.
If your Fs is a scalar, your code should work. But I can't try it out myself right now...
written 7 months ago by klunkean
Thank you very much.
I wrote the code :
print(shape(ps))
print(shape(u_))
​
and I got the answer :
(3,)
(3,)
(3, 2)​
So it appears that fenics doesn't understand that my model is 3D altough my mesh is like this.
written 7 months ago by Liam Gonzalez

3Dmesh

It seems your mesh is surface mesh not volume mesh (please refer to above mesh). Define mesh as volume mesh in GMSH.
written 7 months ago by hirshikesh
the command "print(mesh.geometry().dim())" returns me 2 for this mesh. I do not understand why my final mesh is interpreted as a 2D mesh in my fenics script. look here, it seems to be a volume mesh. Moreover, gmsh tells me that it is meshing volume. Can I force te dimension of the mesh in fenics ?
My .geo script is the following :
d = 0.008 ;
l = 0.4 ;
R = 0.08 ;
e = 0.010;

Point(1) = {0,R,0,d};
Point(2)={0,R+e,0,d};
Point(3) = {l,R+e,0,d};
Point(4) = { l , R , 0 ,d };

Line(1) = { 1,2};
Line(2) = { 2,3};
Line(3) = { 3,4};
Line(4) = { 4,1};
Line loop(1) = {1,2,3,4};
Plane Surface(1) = {1} ;

Extrude {{1, 0, 0}, {0, 0, 0}, Pi/2} {
Surface{1}; Recombine ;
};
Extrude {{1, 0, 0}, {0, 0, 0}, Pi/2} {
Surface{26}; Recombine ;
};
Extrude {{1, 0, 0}, {0, 0, 0}, Pi/2} {
Surface{25}; Surface{48}; Recombine ;
};
Extrude {{1, 0, 0}, {0, 0, 0}, Pi/2} {
Surface{1}; Surface{91}; Recombine ;
};

// definition of the volume
Surface Loop(1) = {39, 43, 35, 47, 25, 90, 78, 82, 86, 104, 108, 112, 100, 17, 21, 13};
Volume(7) = {1};
Physical Volume(1) = {7};

// free surface
Physical Surface(5) = {100, 78, 35, 13, 39, 47, 90, 82, 112, 25};

// submited surface
Physical Surface(6) = {17, 104};

// fixed surface
Physical Surface(7) = {108, 21, 43, 86};
​
And I use dolfin-convert to convert my .msh to a .xml format.
And my first command lines in my fenics script are :
mesh = Mesh("tube.xml")
subdomains = MeshFunction("size_t", mesh, "tube_physical_region.xml")
print(mesh.geometry().dim())​
have a good day.
written 7 months ago by Liam Gonzalez
I finally find out what the problem is.
I use Gmesh3 to produce my mesh. I start with my .geo file and I mesh it (in 3D mode of course) and I visualize this mesh in the graphic interface. However when I re-open the .msh file I got this one.
I really don't understand why. how do you got the mesh you showed me hirshikesh?
written 7 months ago by Liam Gonzalez
Please see the .geo file. print(mesh.geometry().dim()) will give 3
DefineConstant[ lc = { 0.1, Path "Gmsh/Parameters"}];
Point(1) = {0, 0, 0, lc};
Point(2) = {0.5, 0, 0, lc};
Point(3) = {0, 0.5, 0, lc};
Point(4) = {0, -0.5, 0, lc};
Point(5) = {-0.5, 0, 0, lc};
Point(6) = {0.6, 0, 0, lc};
Point(7) = {-0.6, 0, 0, lc};
Point(8) = {0, 0.6, 0, lc};
Point(9) = {0, -0.6, 0, lc};
Circle(1) = {2, 1, 3};
Circle(2) = {3, 1, 5};
Circle(3) = {5, 1, 4};
Circle(4) = {4, 1, 2};
Circle(5) = {9, 1, 6};
Circle(6) = {6, 1, 8};
Circle(7) = {8, 1, 7};
Circle(8) = {7, 1, 9};
Line Loop(9) = {6, 7, 8, 5};
Line Loop(10) = {1, 2, 3, 4};
Plane Surface(11) = {9, 10};
Extrude {0, 0, 1} {
Surface{11};
}
Surface Loop(54) = {48, 11, 24, 28, 32, 36, 53, 40, 44, 52};
Volume(55) = {54};​
written 7 months ago by hirshikesh
Now it's okay, it works well. But I still don't understand the difference beetwin our gmesh scripts, you mix extrusion and direct definition (point, circle and lines) but it works. Nevertheless, thank very much hirshikesh, I'll use your script.
written 7 months ago by Liam Gonzalez

I noticed you are using Identity(2) when you should be using Identity(3). This won't fix the issue with shape(grad(u_)) though. On another note, have you double checked the geometric and topological dimensions of your mesh once it is loaded? E.g.
print(mesh.geometry().dim())
print(mesh.topology().dim())​