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):
    return sym(grad(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
    return sym(grad(w))
  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_))
print(shape(grad(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  
Thank you for all your advise.
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  

1 Answer


0
7 months ago by
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())​
Please login to add an answer/comment or follow this question.

Similar posts:
Search »