Mesh: tradeoff

Return to GFtbox hints and tips

Subdividing the mesh to ensure sufficient resolution to compute curves properly

There is a tradeoff between speed and the accuracy with which curves are computed in meshes. Curves are approximated by a series of straight lines. With too few points, curves will be discontinuous and computed inaccurately. With too many points the computation will be slow. It is, therefore, desirable to vary the density of vertices in the initial mesh according to the expected curvature (you have to run models many times to get this right). This is an example of subdividing a rectangular mesh in the region in which curves will develop.

Simple mesh of 382 elements equally spaced vertices. Superimposed are the polarity arrows (pointing bottom left) and purple and blue factors that control local growth rates.

After one step in which the central region is subdivided twice (the vertices are too dense to see).

As for the previous example but the size of the triangles forming the mesh is only visible along the edge. Red marks regions of increased growth on the top and bottom surfaces.

After growing for 40 steps the excess growth has bent the canvas into an 'S' shape. Note the smooth bend is associated with about 4 elements.

As for the previous example, but there has been no subdivision and the curves are unacceptably incomplete and jagged. There should always be enough elements to permit curves to develop smoothly.

It still might not be correct because the stepsize might be too big, See time tradeoffs It still might not be correct because the tolerances might be too large, See tolerance tradeoffs

Tutorial

Run the project (example project )for 40 steps to see the effect of:

1. subdivision

Then change the modelname to 'NOSUBDIVISION', i.e. make a change to the interaction function

       % Set up names for variant models.  Useful for running multiple models on a cluster.
       m.userdata.ranges.modelname.range = { 'NOSUBDIVISION', 'WITHSUBDIVISION' };  % CLUSTER
       m.userdata.ranges.modelname.index = 2;                       % CLUSTER

by setting index to 1. With too few elements, bends are jagged.

The interaction function is shown below, red highlights the region of interest.

function m = gpt_demosubdivision_20121116( m )
%m = gpt_demosubdivision_20121116( m )
%   Morphogen interaction function.
%   Written at 2012-11-16 12:30:06.
%   GFtbox revision 4351, .

% The user may edit any part of this function between delimiters
% of the form "USER CODE..." and "END OF USER CODE...".  The
% delimiters themselves must not be moved, edited, deleted, or added.

    if isempty(m), return; end

    fprintf( 1, '%s found in %s\n', mfilename(), which(mfilename()) );

    try
        m = local_setproperties( m );
    catch
    end

    realtime = m.globalDynamicProps.currenttime;


%%% USER CODE: INITIALISATION

% In this section you may modify the mesh in any way whatsoever.
    if (Steps(m)==0) && m.globalDynamicProps.doinit % First iteration
        % Zero out a lot of stuff to create a blank slate.  
        % If no morphogens are set in the GUI it may be useful to
        % zero some arrays by uncommenting the following.
        % m.morphogens(:) = 0;
        % m.morphogenclamp(:) = 0;
        % m.mgen_production(:) = 0;
        % m.mgen_absorption(:) = 0;
        % m.seams(:) = false;
        % m.mgen_dilution(:) = false;

        % Set up names for variant models.  Useful for running multiple models on a cluster.
        m.userdata.ranges.modelname.range = { 'NOSUBDIVISION', 'WITHSUBDIVISION' };  % CLUSTER
        m.userdata.ranges.modelname.index = 1;                       % CLUSTER
    end
    modelname = m.userdata.ranges.modelname.range{m.userdata.ranges.modelname.index};  % CLUSTER
    
    disp(sprintf('\nRunning %s model %s\n',mfilename, modelname));    	
    % to plot polariser on the A side and resultant areal growth rate on the B side:
    m = leaf_plotoptions( m, 'morphogenA', 'KAPAR', 'morphogenB', 'KBPAR' ); 
%%% END OF USER CODE: INITIALISATION


%%% SECTION 1: ACCESSING MORPHOGENS AND TIME.
%%% AUTOMATICALLY GENERATED CODE: DO NOT EDIT.

    if isempty(m), return; end

    setGlobals();
    global gNEW_KA_PAR gNEW_KA_PER gNEW_KB_PAR gNEW_KB_PER
    global gNEW_K_NOR gNEW_POLARISER gNEW_STRAINRET gNEW_ARREST
    dt = m.globalProps.timestep;
    polariser_i = gNEW_POLARISER;
    P = m.morphogens(:,polariser_i);
    [kapar_i,kapar_p,kapar_a,kapar_l] = getMgenLevels( m, 'KAPAR' );
    [kaper_i,kaper_p,kaper_a,kaper_l] = getMgenLevels( m, 'KAPER' );
    [kbpar_i,kbpar_p,kbpar_a,kbpar_l] = getMgenLevels( m, 'KBPAR' );
    [kbper_i,kbper_p,kbper_a,kbper_l] = getMgenLevels( m, 'KBPER' );
    [knor_i,knor_p,knor_a,knor_l] = getMgenLevels( m, 'KNOR' );
    [strainret_i,strainret_p,strainret_a,strainret_l] = getMgenLevels( m, 'STRAINRET' );
    [arrest_i,arrest_p,arrest_a,arrest_l] = getMgenLevels( m, 'ARREST' );
    [id_a_i,id_a_p,id_a_a,id_a_l] = getMgenLevels( m, 'ID_A' );
    [id_b_i,id_b_p,id_b_a,id_b_l] = getMgenLevels( m, 'ID_B' );
    [id_subdivide_i,id_subdivide_p,id_subdivide_a,id_subdivide_l] = getMgenLevels( m, 'ID_SUBDIVIDE' );

% Mesh type: rectangle
%            base: 0
%          centre: 0
%      randomness: 0.1
%         version: 1
%           xdivs: 16
%          xwidth: 16
%           ydivs: 8
%          ywidth: 8

%            Morphogen    Diffusion   Decay   Dilution   Mutant
%            --------------------------------------------------
%                KAPAR         ----    ----       ----     ----
%                KAPER         ----    ----       ----     ----
%                KBPAR         ----    ----       ----     ----
%                KBPER         ----    ----       ----     ----
%                 KNOR         ----    ----       ----     ----
%            POLARISER         ----    ----       ----     ----
%            STRAINRET         ----    ----       ----     ----
%               ARREST         ----    ----       ----     ----
%                 ID_A         ----    ----       ----     ----
%                 ID_B         ----    ----       ----     ----
%         ID_SUBDIVIDE         ----    ----       ----     ----


%%% USER CODE: MORPHOGEN INTERACTIONS

% In this section you may modify the mesh in any way that does not

% alter the set of nodes.

    if (Steps(m)==0) && m.globalDynamicProps.doinit  % Initialisation code.
        % Put any code here that should only be performed at the start of
        % the simulation, for example, to set up initial morphogen values.
        id_a_p((m.nodes(:,1)>0)&(m.nodes(:,1)<=2))=1;
        id_b_p((m.nodes(:,1)>-2)&(m.nodes(:,1)<=0))=1;
        id_subdivide_p((m.nodes(:,1)>=-2)&(m.nodes(:,1)<=2))=1;
        P=m.nodes(:,1);
    end
    if realtime<=1
        kapar_p(:) = 0;
        kaper_p(:) = 0;
        kbpar_p(:) = 0;
        kbper_p(:) = 0;
        knor_p(:)  = 0;
    else
        kapar_p(:) = 0.05*id_a_p;
        kaper_p(:) = 0;
        kbpar_p(:) = 0.05*id_b_p;
        kbper_p(:) = 0;
        knor_p(:)  = 0;
    end

%%% END OF USER CODE: MORPHOGEN INTERACTIONS

%%% SECTION 3: INSTALLING MODIFIED VALUES BACK INTO MESH STRUCTURE
%%% AUTOMATICALLY GENERATED CODE: DO NOT EDIT.
    m.morphogens(:,polariser_i) = P;
    m.morphogens(:,kapar_i) = kapar_p;
    m.morphogens(:,kaper_i) = kaper_p;
    m.morphogens(:,kbpar_i) = kbpar_p;
    m.morphogens(:,kbper_i) = kbper_p;
    m.morphogens(:,knor_i) = knor_p;
    m.morphogens(:,strainret_i) = strainret_p;
    m.morphogens(:,arrest_i) = arrest_p;
    m.morphogens(:,id_a_i) = id_a_p;
    m.morphogens(:,id_b_i) = id_b_p;
    m.morphogens(:,id_subdivide_i) = id_subdivide_p;

%%% USER CODE: FINALISATION

% In this section you may modify the mesh in any way whatsoever.
    switch modelname
        case 'NOSUBDIVISION' 
            % do nothing
        case 'WITHSUBDIVISION' 
            % subdivide on the first step
            if realtime>0 && realtime<=0+dt
                m = leaf_subdivide( m, 'morphogen','id_subdivide',...
                    'min',0.5,'max',1,...
                    'mode','mid','levels','all');
            end
            if realtime>0 && realtime<=0+dt
                m = leaf_subdivide( m, 'morphogen','id_subdivide',...
                    'min',0.5,'max',1,...
                    'mode','mid','levels','all');
            end
        otherwise
            % If this happens, maybe you forgot a model.
    end


%%% END OF USER CODE: FINALISATION

end