[Getfem-commits] r5080 - in /trunk/getfem: interface/src/ interface/tests/matlab/ src/

[Yves . Renard]

Author: renard
Date: Fri Sep  4 22:01:14 2015
New Revision: 5080

URL: http://svn.gna.org/viewcvs/getfem?rev=5080&view=rev
Log:
a working example which uses xfem_plus and xfem_minus

Modified:
    trunk/getfem/interface/src/gf_mesh_fem_set.cc
    trunk/getfem/interface/tests/matlab/demo_crack.m
    trunk/getfem/src/getfem_mesh_fem_global_function.cc

Modified: trunk/getfem/interface/src/gf_mesh_fem_set.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/src/gf_mesh_fem_set.cc?rev=5080&r1=5079&r2=5080&view=diff
==============================================================================
--- trunk/getfem/interface/src/gf_mesh_fem_set.cc       (original)
+++ trunk/getfem/interface/src/gf_mesh_fem_set.cc       Fri Sep  4 22:01:14 2015
@@ -25,6 +25,7 @@
 #include <getfem/getfem_partial_mesh_fem.h>
 #include <gmm/gmm_range_basis.h>
 #include <getfem/getfem_mesh_fem_level_set.h>
+#include <getfem/getfem_mesh_fem_product.h>
 
 using namespace getfemint;
 
@@ -259,8 +260,8 @@
        );
 
     /address@hidden ('adapt')
-    For a @tmf levelset object only. Adapt the mesh_fem object to a
-    change of the levelset function. @*/
+      For a @tmf levelset object only. Adapt the mesh_fem object to a
+      change of the levelset function. @*/
     sub_command
       ("adapt", 0, 0, 0, 0,
        getfem::mesh_fem_level_set *mfls
@@ -268,6 +269,19 @@
        if (!mfls) THROW_BADARG("The command 'adapt' can only be "
                                "applied to a mesh_fem_level_set object");
        mfls->adapt();
+       );
+
+    /address@hidden ('set enriched dofs', @ivec DOFs)
+      For a @tmf product object only. Set te enriched dofs and adapt the @tmf 
product.
+      @*/
+    sub_command
+      ("set enriched dofs", 1, 1, 0, 0,
+       getfem::mesh_fem_product *mfprod
+       = dynamic_cast<getfem::mesh_fem_product *>(mf);
+       if (!mfprod) THROW_BADARG("The command 'set enriched dofs' can only be "
+                                 "applied to a mesh_fem_product object");
+       dal::bit_vector doflst = in.pop().to_bit_vector();
+       mfprod->set_enrichment(doflst);
        );
 
   }

Modified: trunk/getfem/interface/tests/matlab/demo_crack.m
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/tests/matlab/demo_crack.m?rev=5080&r1=5079&r2=5080&view=diff
==============================================================================
--- trunk/getfem/interface/tests/matlab/demo_crack.m    (original)
+++ trunk/getfem/interface/tests/matlab/demo_crack.m    Fri Sep  4 22:01:14 2015
@@ -22,11 +22,12 @@
 %
 %%%
 
-variant = 1;
+variant = 4;
 % variant : 1 : one crack with cutoff enrichement
 %           2 : one crack with a fixed size area Xfem enrichment
 %           3 : a branching crack with a fixed size area Xfem enrichment
-%           4 : variant 2 plus a penalisation of the jump over the cracks
+%           4 : variant 3 with the second crack closed by a penalisation of
+%               the jump (exemple of use of xfem_plus and xfem_minus).
 
 
 gf_workspace('clear all');
@@ -67,7 +68,7 @@
 ck3 = gf_global_function('crack',3);
 
 % Definition of the enriched finite element method
-mfls_u    = gf_mesh_fem('levelset',mls,mf_pre_u);
+mfls_u = gf_mesh_fem('levelset',mls,mf_pre_u);
 
 if (variant == 1) % Cutoff enrichement 
   coff = gf_global_function('cutoff',2,0.4,0.01,0.4);
@@ -80,11 +81,29 @@
 else
   mf_part_unity = gf_mesh_fem(m);
   gf_mesh_fem_set(mf_part_unity, 'classical fem', 1);
-  % + selection des ddls autour du ou des fond de fissure
+  DOFpts = gf_mesh_fem_get(mf_part_unity, 'basic dof nodes');
+  % Search the dofs to be enriched with the asymptotic displacement.
+  Idofs_center = find((DOFpts(1,:)).^2 + (DOFpts(2,:)).^2 <= (0.1)^2);
+  mf_sing_u = gf_mesh_fem('global function',m,ls, {ck0,ck1,ck2,ck3}, 1);
+  mf_xfem_sing = gf_mesh_fem('product', mf_part_unity, mf_sing_u);
+  gf_mesh_fem_set( mf_xfem_sing, 'set enriched dofs', Idofs_center);
+  if (variant > 2)
+    Idofs_up = find((DOFpts(1,:)+0.125).^2 + (DOFpts(2,:)-0.375).^2 <= 
(0.1)^2);
+    Idofs_down = find((DOFpts(1,:)+0.125).^2 + (DOFpts(2,:)+0.375).^2 <= 
(0.1)^2);
+    mf_sing_u2 = gf_mesh_fem('global function',m,ls2, {ck0,ck1,ck2,ck3}, 1);
+    mf_xfem_sing2 = gf_mesh_fem('product', mf_part_unity, mf_sing_u2);
+    gf_mesh_fem_set(mf_xfem_sing2, 'set enriched dofs', [Idofs_up Idofs_down]);
+  end
+  
+  if (variant == 2)
+    mf_u = gf_mesh_fem('sum', mf_xfem_sing, mfls_u);
+  else
+    mf_u = gf_mesh_fem('sum', mf_xfem_sing, mf_xfem_sing2, mfls_u);
+  end
 end  
   
 gf_mesh_fem_set(mf_u,'qdim',2);
-% exact solution:
+% Exact solution for a single crack
 mf_ue = gf_mesh_fem('global function',m,ls,{ck0,ck1,ck2,ck3});
 A = 2+2*Mu/(Lambda+2*Mu);
 B=-2*(Lambda+Mu)/(Lambda+2*Mu)
@@ -109,6 +128,12 @@
 
gf_model_set(md,'add_isotropic_linearized_elasticity_brick',mim,'u','lambda','mu');
 gf_model_set(md,'add_initialized_fem_data','DirichletData', mf_ue, Ue);
 gf_model_set(md,'add_Dirichlet_condition_with_penalization',mim,'u', 1e12, 
DIRICHLET, 'DirichletData');
+
+if (variant == 4)
+  mim_bound = gf_mesh_im('levelset', mls, 'boundary(b)', 
gf_integ('IM_STRUCTURED_COMPOSITE(IM_TRIANGLE(6),3)'));
+  % gf_asm('generic', mim_bound, 0, '1', -1)
+  gf_model_set(md, 'add linear generic assembly brick', mim_bound, 
'1e15*(Xfem_plus(u)-Xfem_minus(u)).(Xfem_plus(Test_u)-Xfem_minus(Test_u))');
+end
 
 % assembly of the linear system and solve:
 gf_model_get(md,'solve');

Modified: trunk/getfem/src/getfem_mesh_fem_global_function.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_mesh_fem_global_function.cc?rev=5080&r1=5079&r2=5080&view=diff
==============================================================================
--- trunk/getfem/src/getfem_mesh_fem_global_function.cc (original)
+++ trunk/getfem/src/getfem_mesh_fem_global_function.cc Fri Sep  4 22:01:14 2015
@@ -582,6 +582,8 @@
       update_mls(c.convex_num());
       scalar_type x = mls_x(c.xref());
       scalar_type y = mls_y(c.xref());
+      if (c.xfem_side() > 0 && y <= 0) y = 1E-13;
+      if (c.xfem_side() < 0 && y >= 0) y = -1E-13;
       return fn.val(x,y);
     }
     virtual void grad(const fem_interpolation_context& c,
@@ -591,7 +593,8 @@
       base_small_vector dx(P), dy(P), dfr(2);
       scalar_type x = mls_x.grad(c.xref(), dx);
       scalar_type y = mls_y.grad(c.xref(), dy);
-
+      if (c.xfem_side() > 0 && y <= 0) y = 1E-13;
+      if (c.xfem_side() < 0 && y >= 0) y = -1E-13;
       base_small_vector gfn = fn.grad(x,y);
       gmm::mult(c.B(), gfn[0]*dx + gfn[1]*dy, g);
     }
@@ -603,7 +606,8 @@
       base_small_vector dx(P), dy(P), dfr(2),  dx_real(N), dy_real(N);
       scalar_type x = mls_x.grad(c.xref(), dx);
       scalar_type y = mls_y.grad(c.xref(), dy);
-
+      if (c.xfem_side() > 0 && y <= 0) y = 1E-13;
+      if (c.xfem_side() < 0 && y >= 0) y = -1E-13;
       base_small_vector gfn = fn.grad(x,y);
       base_matrix hfn = fn.hess(x,y);