Dear Yves,
It worked correctly now. Thanks a lot.
Best regards,
Phuoc
On 19/06/2018 15:49, Yves Renard wrote:
Dear Phuoc,
Yes, the export can also have some problems with the singular
transformation ... I robustified a bit the evaluation. Can you
try ?
Best regards,
Yves
Le 19/06/2018 à 12:48, Huu Phuoc BUI a écrit :
Dear Yves,
I now faced with mixed elements of higher order. In the
enclosed script, I use degree = 2 for all elements, and set
the integration method mim = gf.MeshIm(m, 8), but after
resolution, there are 'nan' values in displacement filed. The
'nan' values are shown by yellow color in Paraview as
I tried to change the order for integration method but
without success.
Please let me know if you have an idea. Thanks a lot.
Best regards,
Phuoc
On 19/06/2018 10:21, Huu Phuoc BUI
wrote:
Dear Yves,
Thanks a lot for your help. I think the fact that assigning
the integration method by only providing the degree of
interpolation is very useful. That solved my problem now.
Thank you also for your remark about using macro for
models. It's really helpful.
Best regards,
Phuoc
On 18/06/2018 17:32, Yves Renard
wrote:
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
---------
--
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
---------
|