[Getfem-users] establishing only Neumann and Dirichlet conditions for right side vector(volume forces == 0) in case of xfem dof addition

[Egor Vtorushin]

Hello,  Yves and all

I've been investigating the crack.cc test and trying to set only Neumann and Dirichlet conditions for right side vector. So instead of

<code>
// Defining the volumic source term.
  plain_vector F(nb_dof_rhs * N);
  for (size_type i = 0; i < nb_dof_rhs; ++i)
      gmm::copy(sol_f(mf_rhs.point_of_dof(i)),
        gmm::sub_vector(F, gmm::sub_interval(i*N, N)));
  
  // Volumic source term brick.
  getfem::mdbrick_source_term<> VOL_F(*pINCOMP, mf_rhs, F);
</code>

i write
<code>
// Defining the volumic source term.
  plain_vector F(nb_dof_rhs * N);
      // Volumic source term brick.
  getfem::mdbrick_source_term<> VOL_F(*pINCOMP, mf_rhs, F);
</code>

Thus now the volume forces is set to zero.
Now it is necessary to define track forces on boundary where Neumann condition type is situated. To do this insted of
<code>
 // Defining the Neumann condition right hand side.
  gmm::clear(F);
  
  // Neumann condition brick.
        
  getfem::mdbrick_source_term<>  NEUMANN(VOL_F, mf_rhs, F,NEUMANN_BOUNDARY_NUM);   
</code>
i write
<code>
// Defining the Neumann condition right hand side.
for (size_type i = 0; i < nb_dof_rhs; ++i){
            
        F[i] = 0;
        F[i+1]=0.1*
rhs_func(mf_rhs.point_of_dof(i)[0], mf_rhs.point_of_dof(i)[1]);
}
// Neumann condition brick.

getfem::mdbrick_source_term<> NEUMANN(VOL_F, mf_rhs, F, NEUMANN_BOUNDARY_NUM);
</code>
Here 
rhs_func the auxiliary function equaled to 1 near(to 1e-8) the boundary where Nuemann condition type is situated. But results of program does not corresponded   with BC we just described. So the question - how can i establish the BC i need(only Neumann and Dirichlet conditions for right side vector).
I think the reason that actually number of dof is greater than  vector F size because of xfem addition for describing crack and singularity but i don't now how to  solve this situation. Anybody can help me?  

Regards, Egor Vtorushin

p.s.
the crack.param is given


% -*- matlab -*- (enables emacs matlab mode)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% parameters for program crack                                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%% pde parameters :                                  %%%%%
% MU = 77.0;            % Lamй coefficient.
% LAMBDA = 107.0;       % Lamй coefficient.

MU = 1.0;            % Lamй coefficient.
LAMBDA = 1.0;           % Lamй coefficient.


QUAD = 0;


%%%%%   discretisation parameters  :                               %%%%%

if (~QUAD) 
  MESH_TYPE = 'GT_PK(2,1)';         % linear triangles
else
  % MESH_TYPE = 'GT_LINEAR_QK(2)';
  MESH_TYPE = 'GT_QK(2, 1)';
end;


NX = 15;                          % space step.

MESH_NOISED = 0; % Set to one if you want to "shake" the mesh



if (~QUAD)
  %FEM_TYPE = 'FEM_PK_WITH_CUBIC_BUBBLE(2, 2)';
  FEM_TYPE = 'FEM_PK(2, 1)';  % PK element
  DATA_FEM_TYPE = 'FEM_PK(2,1)';
  INTEGRATION = 'IM_STRUCTURED_COMPOSITE(IM_TRIANGLE(6), 5)';
  %INTEGRATION = 'IM_TRIANGLE(6)';
  FEM_TYPE_P = 'FEM_PK(2,1)';
  MORTAR_FEM_TYPE = FEM_TYPE;
else
  FEM_TYPE = 'FEM_QK(2,2)';  % Q1 fem for quadrangles
  DATA_FEM_TYPE = 'FEM_QK(2,2)';
  INTEGRATION = 'IM_GAUSS_PARALLELEPIPED(2, 5)';
  FEM_TYPE_P = 'FEM_QK(2,1)'; 
  MORTAR_FEM_TYPE = FEM_TYPE;
end;


MIXED_PRESSURE=0;       % Mixed version or not.
DIRICHLET_VERSION = 2;

% integration meth. for sub-simplexe of elements crossed by the level-set
SIMPLEX_INTEGRATION = 'IM_STRUCTURED_COMPOSITE(IM_TRIANGLE(6),3)';

% integration meth. for quasi-polar integration of sub-simplexes adjascent to the level-set
% (comment it to disable quasipolar integration). Should be a
% method defined on a square for 2D, or defined on a prism for 3D.
% SINGULAR_INTEGRATION = 'IM_GAUSS_PARALLELEPIPED(2, 10)';
SINGULAR_INTEGRATION = 'IM_STRUCTURED_COMPOSITE(IM_GAUSS_PARALLELEPIPED(2, 6), 9)';

ADDITIONAL_CRACK = 0;

%Enable the following two lines to use the precalculated solution as enrichement 
%GLOBAL_FUNCTION_MF = "crack.meshfem"
%GLOBAL_FUNCTION_U  = "crack.U"

ENRICHMENT_OPTION = 1;  % 0 = Pas d'enrichissement
                    % 1 = standard XFEM on a fixed zone
            % 2 = global functions with mortar junction
                % 3 = global functions with cutoff
            % 4 = spider fem alone
                    % 5 = spider fem enrichment


RADIUS_ENR_AREA = 0.1;
CUTOFF_FUNC = 2; % 0 for the exponential cutoff. 
                 % 1 for a 3rd degree polynomial cutoff.
                 % 2 for a 5th degree polynomial cutoff. 
CUTOFF = 0.4;
CUTOFF1 = 0.01;
CUTOFF0 = 0.49;

SPIDER_RADIUS =  0.3;
SPIDER_NR = 15;         % size of the cartesian mesh in r for spider fem
SPIDER_NTHETA = 15;      % size of the cartesian mesh in theta for spider fem

SPIDER_K=1;             % order of the  spider fem 


RESIDUAL = 1E-9;         % residual for iterative methods if any.


%%%%%   saving parameters                                             %%%%%
ROOTFILENAME = 'crack';     % Root of data files.
VTK_EXPORT = 1 % export solution to a .vtk file ?