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.

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.

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


%%% USER CODE: SUBFUNCTIONS


function m = local_setproperties( m )
% This function is called at time zero in the INITIALISATION section of the
% interaction function.  It provides commands to set each of the properties
% that are contained in m.globalProps.  Uncomment whichever ones you would
% like to set yourself, and put in whatever value you want.
%
% Some of these properties are for internal use only and should never be
% set by the user.  At some point these will be moved into a different
% component of m, but for the present, just don't change anything unless
% you know what it is you're changing.

%    m = leaf_setproperty( m, 'trinodesvalid', true );
%    m = leaf_setproperty( m, 'prismnodesvalid', true );
%    m = leaf_setproperty( m, 'thicknessRelative', 1.600000 );
%    m = leaf_setproperty( m, 'thicknessArea', 0.000000 );
%    m = leaf_setproperty( m, 'thicknessMode', 'physical' );
%    m = leaf_setproperty( m, 'activeGrowth', 1.000000 );
%    m = leaf_setproperty( m, 'displayedGrowth', 1.000000 );
%    m = leaf_setproperty( m, 'displayedMulti', [] );
%    m = leaf_setproperty( m, 'allowNegativeGrowth', true );
%    m = leaf_setproperty( m, 'usePrevDispAsEstimate', true );
%    m = leaf_setproperty( m, 'perturbInitGrowthEstimate', 0.000010 );
%    m = leaf_setproperty( m, 'perturbRelGrowthEstimate', 0.010000 );
%    m = leaf_setproperty( m, 'perturbDiffusionEstimate', 0.000100 );
%    m = leaf_setproperty( m, 'resetRand', false );
%    m = leaf_setproperty( m, 'mingradient', 0.000000 );
%    m = leaf_setproperty( m, 'relativepolgrad', false );
%    m = leaf_setproperty( m, 'usefrozengradient', true );
%    m = leaf_setproperty( m, 'userpolarisation', false );
%    m = leaf_setproperty( m, 'thresholdsq', 4.000000 );
%    m = leaf_setproperty( m, 'splitmargin', 1.400000 );
%    m = leaf_setproperty( m, 'splitmorphogen',  );
%    m = leaf_setproperty( m, 'thresholdmgen', 0.500000 );
%    m = leaf_setproperty( m, 'bulkmodulus', 1.000000 );
%    m = leaf_setproperty( m, 'unitbulkmodulus', true );
%    m = leaf_setproperty( m, 'poissonsRatio', 0.300000 );
%    m = leaf_setproperty( m, 'starttime', 0.000000 );
%    m = leaf_setproperty( m, 'timestep', 0.010000 );
%    m = leaf_setproperty( m, 'timeunitname',  );
%    m = leaf_setproperty( m, 'distunitname', 'mm' );
%    m = leaf_setproperty( m, 'scalebarvalue', 0.000000 );
%    m = leaf_setproperty( m, 'validateMesh', true );
%    m = leaf_setproperty( m, 'rectifyverticals', false );
%    m = leaf_setproperty( m, 'allowSplitLongFEM', true );
%    m = leaf_setproperty( m, 'longSplitThresholdPower', 0.000000 );
%    m = leaf_setproperty( m, 'allowSplitBentFEM', false );
%    m = leaf_setproperty( m, 'allowSplitBio', true );
%    m = leaf_setproperty( m, 'allowFlipEdges', false );
%    m = leaf_setproperty( m, 'allowElideEdges', true );
%    m = leaf_setproperty( m, 'mincellangle', 0.200000 );
%    m = leaf_setproperty( m, 'alwaysFlat', 0.000000 );
%    m = leaf_setproperty( m, 'flattenforceconvex', true );
%    m = leaf_setproperty( m, 'flatten', false );
%    m = leaf_setproperty( m, 'flattenratio', 1.000000 );
%    m = leaf_setproperty( m, 'useGrowthTensors', false );
%    m = leaf_setproperty( m, 'plasticGrowth', false );
%    m = leaf_setproperty( m, 'totalinternalrotation', 0.000000 );
%    m = leaf_setproperty( m, 'stepinternalrotation', 2.000000 );
%    m = leaf_setproperty( m, 'showinternalrotation', false );
%    m = leaf_setproperty( m, 'performinternalrotation', false );
%    m = leaf_setproperty( m, 'internallyrotated', false );
%    m = leaf_setproperty( m, 'maxFEcells', 0 );
%    m = leaf_setproperty( m, 'inittotalcells', 0 );
%    m = leaf_setproperty( m, 'bioApresplitproc',  );
%    m = leaf_setproperty( m, 'bioApostsplitproc',  );
%    m = leaf_setproperty( m, 'maxBioAcells', 0 );
%    m = leaf_setproperty( m, 'colors', (6 values) );
%    m = leaf_setproperty( m, 'colorvariation', 0.050000 );
%    m = leaf_setproperty( m, 'colorparams', (12 values) );
%    m = leaf_setproperty( m, 'biocolormode', 'auto' );
%    m = leaf_setproperty( m, 'freezing', 0.000000 );
%    m = leaf_setproperty( m, 'canceldrift', false );
%    m = leaf_setproperty( m, 'mgen_interaction',  );
%    m = leaf_setproperty( m, 'mgen_interactionName', 'gpt_demosubdivision_20121116' );
%    m = leaf_setproperty( m, 'allowInteraction', true );
%    m = leaf_setproperty( m, 'interactionValid', true );
%    m = leaf_setproperty( m, 'gaussInfo', (unknown type struct) );
%    m = leaf_setproperty( m, 'stitchDFs', [] );
%    m = leaf_setproperty( m, 'D', (36 values) );
%    m = leaf_setproperty( m, 'C', (36 values) );
%    m = leaf_setproperty( m, 'G', (6 values) );
%    m = leaf_setproperty( m, 'solver', 'cgs' );
%    m = leaf_setproperty( m, 'solverprecision', 'double' );
%    m = leaf_setproperty( m, 'solvertolerance', 0.001000 );
%    m = leaf_setproperty( m, 'solvertolerancemethod', 'max' );
%    m = leaf_setproperty( m, 'diffusiontolerance', 0.000010 );
%    m = leaf_setproperty( m, 'allowsparse', true );
%    m = leaf_setproperty( m, 'maxIters', 0 );
%    m = leaf_setproperty( m, 'maxsolvetime', 1000.000000 );
%    m = leaf_setproperty( m, 'cgiters', 0 );
%    m = leaf_setproperty( m, 'simsteps', 0 );
%    m = leaf_setproperty( m, 'stepsperrender', 0 );
%    m = leaf_setproperty( m, 'growthEnabled', true );
%    m = leaf_setproperty( m, 'diffusionEnabled', true );
%    m = leaf_setproperty( m, 'flashmovie', false );
%    m = leaf_setproperty( m, 'makemovie', false );
%    m = leaf_setproperty( m, 'moviefile',  );
%    m = leaf_setproperty( m, 'codec', 'None' );
%    m = leaf_setproperty( m, 'autonamemovie', true );
%    m = leaf_setproperty( m, 'overwritemovie', false );
%    m = leaf_setproperty( m, 'framesize', [] );
%    m = leaf_setproperty( m, 'mov', [] );
%    m = leaf_setproperty( m, 'boingNeeded', false );
%    m = leaf_setproperty( m, 'initialArea', 128.000000 );
%    m = leaf_setproperty( m, 'bendunitlength', 11.313708 );
%    m = leaf_setproperty( m, 'targetRelArea', 1.000000 );
%    m = leaf_setproperty( m, 'defaultinterp', 'min' );
%    m = leaf_setproperty( m, 'readonly', false );
%    m = leaf_setproperty( m, 'projectdir', 'D:\ab\Matlab stuff\Growth models' );
%    m = leaf_setproperty( m, 'modelname', 'GPT_DemoSubdivision_20121116' );
%    m = leaf_setproperty( m, 'allowsave', true );
%    m = leaf_setproperty( m, 'addedToPath', false );
%    m = leaf_setproperty( m, 'bendsplit', 0.300000 );
%    m = leaf_setproperty( m, 'usepolfreezebc', false );
%    m = leaf_setproperty( m, 'dorsaltop', true );
%    m = leaf_setproperty( m, 'defaultazimuth', -45.000000 );
%    m = leaf_setproperty( m, 'defaultelevation', 33.750000 );
%    m = leaf_setproperty( m, 'defaultroll', 0.000000 );
%    m = leaf_setproperty( m, 'defaultViewParams', (unknown type struct) );
%    m = leaf_setproperty( m, 'comment',  );
%    m = leaf_setproperty( m, 'legendTemplate', '%T: %q\n%m' );
%    m = leaf_setproperty( m, 'bioAsplitcells', true );
%    m = leaf_setproperty( m, 'bioApullin', 0.142857 );
%    m = leaf_setproperty( m, 'bioAfakepull', 0.202073 );
%    m = leaf_setproperty( m, 'interactive', false );
%    m = leaf_setproperty( m, 'coderevision', 0 );
%    m = leaf_setproperty( m, 'coderevisiondate',  );
%    m = leaf_setproperty( m, 'modelrevision', 0 );
%    m = leaf_setproperty( m, 'modelrevisiondate',  );
%    m = leaf_setproperty( m, 'savedrunname',  );
%    m = leaf_setproperty( m, 'savedrundesc',  );
%    m = leaf_setproperty( m, 'vxgrad', (108 values) );
%    m = leaf_setproperty( m, 'lengthscale', 16.000000 );
end


% Here you may write any functions of your own, that you want to call from
% the interaction function, but never need to call from outside it.
% Remember that they do not have access to any variables except those
% that you pass as parameters, and cannot change anything except by
% returning new values as results.
% Whichever section they are called from, they must respect the same
% restrictions on what modifications they are allowed to make to the mesh.

% For example:


% function m = do_something( m )
%   % Change m in some way.
% end

% Call it from the main body of the interaction function like this:
%       m = do_something( m );