You are an expert in writing code in many different languages. Your goal is to perform code completion for the following code, keeping in mind the rest of the code and the file meta data.
    Metadata details:
    filename: ex_nonnewtonian1.m
    path of file: ./repos/featool-multiphysics/examples
    
    the code of the file until where you have to start completion: function [ fea, out ] = ex_nonnewtonian1( varargin )
%EX_NONNEWTONIAN1 2D Example for non-Newtonian flow in a channel.
%
%   [ FEA, OUT ] = EX_NONNEWTONIAN1( VARARGIN ) Sets up and solves
%   stationary non-Newtonian flow in a rectangular channel. The fluid
%   is governed by a power-law model such that the viscosity can be
%   expressed as mu_eff = mu0*K*D^(n-1). A constant pressure drop
%   is prescribed resulting in an analytical solution for the velocity
%   u(y) = n/(n+1)*(1/(mu0*K)*dp)^(1/n)*((h/2)^((n+1)/n)-abs(h/2-y)^((n+1)/n)).
%   Accepts the following property/value pairs.
%
%       Input       Value/{Default}        Description

    
function [ fea, out ] = ex_nonnewtonian1( varargin )
%EX_NONNEWTONIAN1 2D Example for non-Newtonian flow in a channel.
%
%   [ FEA, OUT ] = EX_NONNEWTONIAN1( VARARGIN ) Sets up and solves
%   stationary non-Newtonian flow in a rectangular channel. The fluid
%   is governed by a power-law model such that the viscosity can be
%   expressed as mu_eff = mu0*K*D^(n-1). A constant pressure drop
%   is prescribed resulting in an analytical solution for the velocity
%   u(y) = n/(n+1)*(1/(mu0*K)*dp)^(1/n)*((h/2)^((n+1)/n)-abs(h/2-y)^((n+1)/n)).
%   Accepts the following property/value pairs.
%
%       Input       Value/{Default}        Description
%       -----------------------------------------------------------------------------------
%       rho         scalar {1}             Density
%       mu0         scalar {1.29684}       Newtonian viscosity
%       dp          scalar {1}             Pressure drop
%       n           scalar {0.25}          Power-law exponent
%       K           scalar {1}             Power-law constant
%       h           scalar {1}             Channel height
%       l           scalar {0.2}           Channel length
%       igrid       scalar 1/{0}           Cell type (0=quadrilaterals, 1=triangles)
%       hmax        scalar {0.1}           Max grid cell size
%       sf_u        string {sflag1}        Shape function for velocity
%       sf_p        string {sflag1}        Shape function for pressure
%       iplot       scalar 0/{1}           Plot solution and error (=1)
%                                                                                         .
%       Output      Value/(Size)           Description
%       -----------------------------------------------------------------------------------
%       fea         struct                 Problem definition struct
%       out         struct                 Output struct

% Copyright 2013-2023 Precise Simulation, Ltd.


cOptDef = { ...
  'rho',      1;
  'mu0',      1.29684;
  'K',        1;
  'n' ,       0.25;
  'dp' ,      1;
  'h',        1;
  'l',        2;
  'igrid',    1;
  'hmax',     0.05;
  'sf_u',     'sflag2';
  'sf_p',     'sflag1';
  'iplot',    1;
  'fid',      1 };
[got,opt] = parseopt(cOptDef,varargin{:});
fid       = opt.fid;


% Geometry and grid parameters.
h         = opt.h;       % Height of rectangular domain.
l         = opt.l;       % Length of rectangular domain.
% Discretization parameters.
sf_u      = opt.sf_u;    % FEM shape function type for velocity.
sf_p      = opt.sf_p;    % FEM shape function type for pressure.


% Geometry definition.
gobj = gobj_rectangle( 0, l, 0, h );
fea.geom.objects = { gobj };
fea.sdim = { 'x' 'y' };   % Coordinate names.


% Grid generation.
if ( opt.igrid==1 )
  fea.grid = gridgen(fea,'hmax',opt.hmax,'fid',fid);
else
  fea.grid = rectgrid(round(l/opt.hmax),round(h/opt.hmax),[0 l;0 h]);
  if( opt.igrid<0 )
    fea.grid = quad2tri( fea.grid );
  end