Re: [Getfem-users] Elastic interface

[Yves Renard]

Dear Jean-François,

Yes, the better is to use simple mesh_fems in that case. The mesh_fem_level_set 
object directly build a field which is discontinuous through the interface. 
But, the problem is that you will have no mean to access to the jump between 
the displacements on the two sides of the interface inside the generic assembly.

Concerning the unit normal vector, it is indeed not defined on the interface.  
But you can easily define it by adding the level-set field as a data of the 
model and using "Grad_lambda/Norm(Grad_lambda)" as a unit normal vector if 
"lambda" is level-set field.

Best regards,

Yves.

----- Original Message -----
From: "Jean-François Barthélémy" <address@hidden>
To: "Yves Renard" <address@hidden>
Cc: address@hidden
Sent: Tuesday, December 30, 2014 11:33:19 PM
Subject: Re: [Getfem-users] Elastic interface

Dear Yves,

Thank you very much for your explanations which helped me to solve my
problem, at least using the first method.

Indeed it works fine with the method based on a single mesh with two
regions such that the interface is defined by convex edges.

Nevertheless I still have some problems with the "levelset" method.

First, I did not used the mesh_fem_level_set (or MeshFem('levelset',...) in
Python) but a simple meshfem on the original regular mesh and partial
meshfems defined on the dofs given by the mesh_im_level_set (outside and
inside the levelset thanks to dof_from_im) to build two differents fields.
It seems that the same is done in demo_fictitious_domain. I did not succeed
in exploiting the mesh_fem_level_set in this context of two fields. Is it a
problem or did you mean something else ?

Secondly, the major problem is that I can't really build the interface
elasticity stiffness with different tangential and normal stiffness because
the add_linear_generic_assembly_brick can't read the "Normal" keyword
(probably because the integration over edges is defined on a levelset and
not on a "real" boundary ?).
The error message is
************
RuntimeError: (Getfem::InterfaceError) -- Error in
getfem_generic_assembly.cc, line 2314 :
Invalid outward unit normal vector. Possible reasons: not on boundary or
transformation failed.
************
To explain more in details, I build the interface integration method by :
mim_bound = gf.MeshIm('levelset',mls,'boundary', gf.Integ('IM_TRIANGLE(3)'))
Then I call the brick :
md.add_linear_generic_assembly_brick(mim_bound,'0.5*(umat-uinc).((kn*(address@hidden
)+kt*(Id(qdim(umat))address@hidden))')
and this boils down to the error message hereabove.
The program works fine if I remove any reference to "Normal", writing for
example '0.5*(umat-uinc).(umat-uinc)' (ie K=Id), or if I manually
reconstruct the Normal for a simple inclusion ( (X-Center)/Norm(X-Center)
for a disk).

Is there a way to circumvent this problem ?

Thank you very much and happy new year in advance.

Best regards,
Jean-François Barthélémy


2014-12-29 20:26 GMT+01:00 Yves Renard <address@hidden>:

>
>
> Dear Jean-François,
>
> The "interpolate tranformation" tool is a generic tool but it can indeed
> be time consuming because for each Gauss point on the boundary, a search of
> the corresponding element on the other mesh is done. Although this search
> is optimized (rtree structures), it can be of course far more expensive
> than a standard assembly (however, in 2D it is done on a 1D interface so
> that it should be relatively fast).
>
> However, in your case, the simpler way should be to do the job with a
> single mesh defining two different regions. You can define two different
> fields on the two regions by restricting the same "mesh_fem" to the dof of
> the corresponding region (with the partial_mesh_fem object, for instance,
> or with a filtered variable of the model object).
> Then, you have to define two integration methods, one on each region to
> add the two different elasticity terms. Then, you can use the generic
> assembly without "interpolate transformation" to add your elasticity term
> on the interface (you have to define a mesh region by selecting the faces
> corresponding to your interface on only one side of the interface).
>
> Another possibility is indeed to use a level-set and the
> mesh_fem_level_set object which describe a finite element space cut by the
> level-set. However, there is no tool for the moment to deal directly with
> the jump at the interface in the generic assembly (I agree that this would
> be interesting in many situations). So that, the easier way is again two
> define two different fields and use the level-set tools of getfem to define
> cut integrations methods (mesh_im_level_set) still using the generic
> assembly to add your elasticity term on the interface (using a integration
> method on the level-set, also provided by mesh_im_level_set). All is
> available on the python/matlab/scilab interface. You have some examples in
> the interface test programs:
>
> interface/tests/python/demo_fictitious_domains.py
> interface/tests/matlab/demo_fictitious_domains.m
> interface/tests/matlab/demo_fictitious_domains_laplacian.m
> interface/tests/matlab/demo_structural_optimization.m
>
> Best regards,
>
> Yves.
>
>
> ----- Original Message -----
> From: "Jean-François Barthélémy" <address@hidden>
> To: address@hidden
> Sent: Sunday, December 28, 2014 12:15:47 PM
> Subject: [Getfem-users] Elastic interface
>
> Dear Getfem users,
>
> I am trying to find an easy and efficient way to take into account an
> elastic linear interface between two different materials in Getfem (using
> the python interface and/or in C++). The problem is defined on a uniform
> elastic matrix containing an elastic inclusion (with different moduli from
> that of the matrix) such that the matrix/inclusion interface is ruled by a
> relationship of the form T=K.[u] where T is the stress vector acting on the
> interface, K is the interface stiffness and [u] is the displacement jump
> (simple bilateral contact). Some finite element codes address this problem
> by allowing to build "joint elements".
>
> I have already succeeded in solving this problem but in a very inefficient
> way I think. I have indeed defined two differents meshes, one for the
> matrix and one for the inclusion. Although occupying the same geometrical
> domain, the interface boundaries are different from one mesh to the other
> (pairs of points at the same place but one in the matric mesh and the other
> one in the inclusion). I then defined two MeshFems, two variables (umat and
> uinc), two integration methods,... and the contribution of the interface to
> the global elastic stiffness "(umat-uinc).K.(Test_umat-Test_uinc)" by using
> the "interpolate transformation" to project fields expressed on one mesh to
> the other one (kind of mortar method on consistent boundaries). It works
> fine on a small mesh in 2D but is too time consuming on large meshes or in
> 3D.
>
> Is there a better way to do it either by using one or two meshes in which
> the interface correspond to element edges or, even better, by resorting to
> a levelset to define the interface independently of the mesh and xfem ?
>
> In the case of a levelset, how is it possible to build the discontinuous
> field of elastic moduli (the elements crossed by the levelset contain two
> different materials but how to correctly define all the degrees of freedom
> including those related to the Heaviside function) ? And finally how can I
> define the interface contribution to the stiffness ie integral_{interface}
> [ (umat-uinc).K.(Test_umat-Test_uinc) ] dS, possibly in a high-level
> generic assembly procedure in Python since I don't know how to access to
> discontinuity terms in the high-level language ?
>
> I haven't really found solutions to my questions in the examples (either
> related to mortar, cracks, contacts...) but I may have missed something.
>
> Thank you very much in advance for any help or piece of code.
>
> Best regards,
> Jean-François Barthélémy
>
>
> _______________________________________________
> Getfem-users mailing list
> address@hidden
> https://mail.gna.org/listinfo/getfem-users
>