两相障碍问题的FEM近似及其后验误差估计(Matlab代码实现）

 💥💥💞💞欢迎来到本博客❤️❤️💥💥

🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。

⛳️座右铭：行百里者，半于九十。

目录

💥1 概述

📚2 运行结果

🎉3 参考文献

🌈4 Matlab代码及详细文章

💥1 概述

文献来源：

对于两相障碍物问题，我们推导出基本误差恒等式，该恒等式产生了到精确解的距离的自然度量。对于这个度量，我们推导出一个可计算的多数函数，对允许的（能量）函数类中的任何函数都有效。证明当且仅当函数与最小化器重合时，多数消失。结果表明，相应的估计没有间隙，因此可以以任何所需的精度评估任何近似的精度。 

📚2 运行结果

 

 

 

 

 部分代码：

add_paths
define_global_variables

testing_on=0;   %more figures than in the paper
labels_on=1;    %figures labels on/off
elements_frames_on=0;           
dirichlet_nodes_show=0;   
tricont_level_zero_value=0.0001;    %level set value for the free boundary

%%%%%%%%% START OF SETUP %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
example=1;                                     % 1 or 2 implemented

iterations_majorant=1e1; 
iterations_majorant_to_printout=[1:10 20:10:100 200:100:iterations_majorant];

levels_energy_error=5;                         %maximal level for the diference of energies
levels_majorant=levels_energy_error;           %maximal level for majorant computation 
level_to_visualize=levels_energy_error;        %level to visualize 

if example==1
    level_exact_energy=levels_energy_error;        %level to estimate the exact energy for example 2 (known for example 1!)
else
    level_exact_energy=levels_energy_error+1;        %level to estimate the exact energy for example 2 (known for example 1!)
end
  
testfolder = 'test2D';     % specify the test folder
cd(testfolder)
    %load mesh;                 % load initial mesh
    f = @loadfun;              % right hand side
%    u = @exact;                % exact solution
cd ..

if example==1   %Example 1D computed using 2D    
    alpha_plus=8;   alpha_minus=8;
    
    nodes2coord = [0 0; 1 0; 2 0; 2 1; 1 1; 0 1];
    nodes2coord = nodes2coord - kron(ones(size(nodes2coord,1),1),[1 0]);
    elems2nodes = [1 2 5; 5 6 1; 2 3 4; 2 4 5];
    dirichlet   = [3 4; 6 1];
    neumann=[1 2; 2 3; 4 5; 5 6;];
    
    if 0
        nodes2coord = [0 0; 2 0; 2 1; 0 1; 1 1/2];
        nodes2coord = nodes2coord - kron(ones(size(nodes2coord,1),1),[1 0]);
        elems2nodes = [1 2 5; 2 3 5; 3 4 5; 4 1 5];
        dirichlet   = [1 2; 2 3; 4 1];
        neumann=[];
    end
       
    C_Omega=2/pi;  %Friedrich's constant
    
    energy_exact=5+1/3;
else %Example F. Bozorgnia      
    alpha_plus=4; alpha_minus=4;
        
    nodes2coord = [0 0; 2 0; 2 2; 0 2; 1 1];
    nodes2coord = nodes2coord - kron(ones(size(nodes2coord,1),1),[1 1]);
    elems2nodes = [1 2 5; 2 3 5; 3 4 5; 4 1 5];
    dirichlet   = [1 2; 2 3; 3 4; 4 1];
    neumann=[];
    
    energy_exact=13;  %this value is only expected, maybe wrong
    
    C_Omega=sqrt(2)/pi;    %Friedrich's constant
end

%%%%%%%%% END OF SETUP %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%finding box parameters
x_min=min(nodes2coord(:,1)); x_max=max(nodes2coord(:,1));
y_min=min(nodes2coord(:,2)); y_max=max(nodes2coord(:,2));

for level=0:level_exact_energy
    fprintf('  Level:  %d \n', level)    
    
    % uniform refinement
    if (level>0)
        tic     
        [nodes2coord,elems2nodes,dirichlet,neumann] = refinement_uniform_2D(nodes2coord,elems2nodes,dirichlet,neumann);
        time_ref=toc;
        fprintf('  Mesh refined: [%f]\n',time_ref)     
    end
    
    %calculate affine transformations for elements
    tic;
    [B_K,b_K,B_K_det] = affine_transformations(nodes2coord,elems2nodes);
    [elems2edges, edges2nodes] = get_edges(elems2nodes);  
    signs             = signs_edges(elems2nodes);
    time_at = toc;
    fprintf('  Calculating affine transf., edge data, and edge signs: [%f]\n',time_at)

    nelems = size(elems2nodes,1);
    nnodes = size(nodes2coord,1);
    nedges = max(max(elems2edges));
    
    all_number_of_nodes(level+1)=size(nodes2coord,1);
    all_number_of_elements(level+1)=size(elems2nodes,1);

    fprintf('  NOF elements = %d, NOF nodes =  %d, NOF edges = %d.\n', nelems, nnodes, nedges)

    % stiffness matrix assembly P1
    tic; K = stiffness_matrix_P1(elems2nodes,B_K,B_K_det); time1 = toc;
    fprintf('  Assembling P1 stifness matrix K [%f], ',time1);

    % mass matrix assembly
    tic; M = mass_matrix_P1(elems2nodes,B_K_det); time2 = toc;
    fprintf('P1 mass matrix M [%f]\n',time2);

    %extracts nodes from Dirichlet BC
    nodes_dirichlet=unique(dirichlet);
    nodes_dirichlet_left=nodes_dirichlet(nodes2coord(nodes_dirichlet,1)==x_min);
    nodes_dirichlet_right=nodes_dirichlet(nodes2coord(nodes_dirichlet,1)==x_max);
    nodes_dirichlet_bottom=nodes_dirichlet(nodes2coord(nodes_dirichlet,2)==y_min);
    nodes_dirichlet_top=nodes_dirichlet(nodes2coord(nodes_dirichlet,2)==y_max);

    %defining Dirichlet and free nodes
    if example==1
        nodes_dirichlet=unique([nodes_dirichlet_left; nodes_dirichlet_right]);     
        v_dirichlet=zeros(nnodes,1);
        v_dirichlet(nodes_dirichlet_left,:)=-ones(size(nodes_dirichlet_left,1),1);
        v_dirichlet(nodes_dirichlet_right,:)=ones(size(nodes_dirichlet_right,1),1);
        %v_dirichlet(nodes_dirichlet_top,:)=(nodes2coord(nodes_dirichlet_top,1)).^3;
        %v_dirichlet(nodes_dirichlet_bottom,:)=(nodes2coord(nodes_dirichlet_bottom,1)).^3;
    else
        v_dirichlet=zeros(nnodes,1);
        v_dirichlet(nodes_dirichlet_left,:)= nodes2coord(nodes_dirichlet_left,2)-1;
        v_dirichlet(nodes_dirichlet_right,:)=nodes2coord(nodes_dirichlet_right,2)+1;
        v_dirichlet(nodes_dirichlet_top,:)=nodes2coord(nodes_dirichlet_top,1)+1;
        v_dirichlet(nodes_dirichlet_bottom,:)=nodes2coord(nodes_dirichlet_bottom,1)-1;
    end
    nodes_free=setdiff((1:nnodes)',nodes_dirichlet);
    
    %define Neumann edges
    edges_neumann=get_boundary_edges(edges2nodes,neumann);
    
    % extracts matrices from K
    K_ii=K(nodes_free,nodes_free); 
    K_dd=K(nodes_dirichlet,nodes_dirichlet);
    K_id=K(nodes_free,nodes_dirichlet);
    
    if level==level_to_visualize
        figure(1);
        show_mesh(elems2nodes,nodes2coord); 
        if labels_on
            if dirichlet_nodes_show
               hold on                    
               plot(nodes2coord(nodes_dirichlet,1),nodes2coord(nodes_dirichlet,2),'.r', 'MarkerSize',15)                
               hold off              
                title('visualization mesh + Dirichlet nodes (red circles)');
            else
                title('visualization mesh');
            end                
            view(2);
        end
        axis equal; axis tight;
    end

🎉3 参考文献

部分理论来源于网络，如有侵权请联系删除。

1]Repin, S., Valdman, J. A Posteriori Error Estimates for Two-Phase Obstacle Problem. J Math Sci 207, 324–335 (2015). https://doi.org/10.1007/s10958-015-2374-9

🌈4 Matlab代码及详细文章