Dear Yves, Thank you very much. It was indeed the most efficient solution and it works fine.
I wish you a happy new year and thanks again for this very nice tool. Best regards, Jean-François Barthélémy 2014-12-31 10:55 GMT+01:00 Yves Renard <[email protected]>: > > > 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" <[email protected]> > To: "Yves Renard" <[email protected]> > Cc: [email protected] > 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*(Normal@Normal > )+kt*(Id(qdim(umat))-Normal@Normal)).(umat-uinc))') > 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 <[email protected]>: > > > > > > > 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" <[email protected]> > > To: [email protected] > > 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 > > [email protected] > > https://mail.gna.org/listinfo/getfem-users > > >
_______________________________________________ Getfem-users mailing list [email protected] https://mail.gna.org/listinfo/getfem-users
