Dear Phuoc,
The main reason is that you use a Newton-Cotes integration methods
on the prism. The Newton-Cotes integration method has Gauss points
on the geometrical nodes which is not convenient at all for the
pyramid. You can use instead
mim = gf.MeshIm(m, 3)
which will select an integration of order at least 3 on each
element.
Just an additional remark. You can use the macro system which can
give clearer expressions. For instance
sig_u = "(lambda_para*Trace(Grad_u)*Id(qdim(u)) + mu_para*(Grad_u
+ Grad_u'))"
stress_jump_inner = "{sig_u} .... ".format(sig_u=sig_u)
can also be entered has
md.add_macro('sig_u', "(lambda_para*Trace(Grad_u)*Id(qdim(u)) +
mu_para*(Grad_u + Grad_u'))")
stress_jump_inner = "sig_u ... "
or
stress_jump_inner = "Def sig_u :=
(lambda_para*Trace(Grad_u)*Id(qdim(u)) + mu_para*(Grad_u +
Grad_u')); sig_u ..."
you can also define macros with parameters, see
http://getfem.org/userdoc/gasm_high.html#macro-definition
Best regards,
Yves
Le 18/06/2018 à 13:35, Huu Phuoc BUI a écrit :
Dear Yves,
Please find enclosed the smallest python script. I made an
mistake when assigning Lamé coefficients for 3D problems
(sorry!). Now the problem seems to be different than the 'nan'
ouput.
By using generic assembly, it says:
"RuntimeError: (Getfem::InterfaceError) -- Error in
bgeot_geometric_trans.cc, line 245 bgeot::scalar_type
bgeot::lu_inverse(bgeot::scalar_type*, bgeot::size_type, bool):
Non invertible matrix
"
Best regards,
Phuoc
On 18/06/2018 10:57, Yves Renard
wrote:
Dear Phuoc,
Do you compute some gradient quantities on the finite element
nodes ? If not, it should not pose any problem since using the
generic assembly, the computation are located on the Gauss
points which are not located on the nodes.
can
you extract the smallest possible program on which the
error occurs ?
Best regards,
Yves
Le 18/06/2018 à 10:35, Huu Phuoc BUI a écrit :
Dear
Yves,
Sorry I email you again on this topics. In fact, the generic
assembly did some computation on pyramid elements but the
output is 'nan'. On elements of others types (of the mixed
mesh), the output is fine.
Please let me know if you have an idea. Thanks.
Best regards,
Phuoc
On 15/06/2018 16:43, Huu Phuoc BUI wrote:
Dear Yves,
That's great! It's indeed working now. Thank you very much
for your help.
Best regards,
Phuoc
On 15/06/2018 16:38, Yves Renard wrote:
Dear Phuoc,
The problem is indeed the same, to compute the location of
gauss points on a neighbour element, the inversion of the
geometric transformation is called. This inversion used a
Newton algorithm whose starting point is a node of the
element ... (function find_initial_guess in
bgeot_geotrans_inv.cc).
I changed a bit the starting point, so it should work now.
But this is clear that the fact that the transformation is
singular at the pyramid tip is a difficulty.
Best regards,
Yves
Le 15/06/2018 à 11:53, Huu Phuoc BUI a écrit :
Dear Yves,
Thank you very much for your help.
Now it worked when computing the outer faces. But I
faced another problem of this mixed mesh. In fact, I'd
like to compute the stress jump at the inner faces using
generic assembly, and GetFEM++ showed an error like
this:
"
File
"/usr/local/lib/python2.7/dist-packages/getfem/getfem.py",
line 5591, in asm_generic
return getfem('asm', 'generic', mim, order,
_expression_, region, model, *args)
RuntimeError: (Getfem::InterfaceError) -- Error in
getfem_generic_assembly_compile_and_exec.cc, line 3887
virtual int
getfem::ga_instruction_neighbour_transformation_call::exec():
Geometric transformation inversion has failed in
neighbour transformation
"
So, as you said, this error comes from geometric
transformation of pyramid elements. Do you have an idea
how to fix it?
Best regards,
Phuoc
On 15/06/2018 11:10, Yves Renard wrote:
Dear Phuoc,
The problem come from the fact that the test is made
on the mean unit normal vector on a face. This mean
unit normal vector is computed by averaging the unit
normal vector on each geometrical node of a face.
Unfortunately, for a pyramid element, the
transformation is singular on the geometrical nodes.
I replaced the mean of normal vector by the normal
vector on the mean of geometrical nodes of the face.
It fixes the problem. (I push the change on the
master branch).
Best regards,
Yves
Le 14/06/2018 à 15:40, Huu Phuoc BUI a écrit :
Dear Yves,
Please find the mesh attached.
Thank you for your help.
Best regards,
Phuoc
On 14/06/2018 14:52, Yves Renard wrote:
Dear Phuoc,
No, I see no reason. If you can just send me the
mesh in a format readable by GetFEM, I could check
what is wrong.
Best regards,
Yves
Le 14/06/2018 à 10:24, Huu Phuoc BUI a écrit :
Dear Getfem-users,
I have a problem with computing the outer faces
when using a mixed mesh (having tetras, hexa,
prisms and pyramids).
The error is:
face_left = m.outer_faces_with_direction([0, 0.,
-1.0], 0.01)
File
"/usr/local/lib/python2.7/dist-packages/getfem/getfem.py",
line 1372, in outer_faces_with_direction
return
self.get("outer_faces_with_direction", v, angle,
CVIDs)
File
"/usr/local/lib/python2.7/dist-packages/getfem/getfem.py",
line 1152, in get
return getfem('mesh_get',self.id, *args)
RuntimeError: (Getfem::InterfaceError) -- Error
in bgeot_geometric_trans.cc, line 381 const
base_matrix&
bgeot::geotrans_interpolation_context::B()
const:
Non invertible matrix
Please let me know if you have any idea. Thank
you.
Best regards,
Phuoc
--
Yves Renard (address@hidden) tel : (33) 04.72.43.87.08
Pole de Mathematiques, INSA-Lyon fax : (33) 04.72.43.85.29
20, rue Albert Einstein
69621 Villeurbanne Cedex, FRANCE
http://math.univ-lyon1.fr/~renard
---------
--
Yves Renard (address@hidden) tel : (33) 04.72.43.87.08
Pole de Mathematiques, INSA-Lyon fax : (33) 04.72.43.85.29
20, rue Albert Einstein
69621 Villeurbanne Cedex, FRANCE
http://math.univ-lyon1.fr/~renard
---------
|