title: How to import non-Newtonian model in FEATOOL Multiphysics?
body: <p>I am using FEATool Multiphysics. I have a MATLAB function named ex_nonnewtonian1.m, which is supported by FEATool itself and contains the definition of a power-law non-Newtonial model. This file wasn't added to FEATOOL, so I need to add it manually.</p>
  
  <p>How can I add the function <code>ex_nonnewtonian1</code> to  the equation box in FEATool window?</p>
  
  <p>This is the code for <code>ex_nonnewtonian1</code>:</p>
  
  <pre><code>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-2020 Precise Simulation, Ltd.
  
  
  cOptDef = { ...   'rho',      1;   'mu0',      1.29684;   'dp' ,      1;   'n' ,       0.25;   'K',        1;   'h',        1;   'l',        2;   'igrid',    1;   'hmax',     0.1;   'sf_u',     'sflag2';   'sf_p',     'sflag1';   'iplot',    1;   'fid',      1 }; [got,opt] = parseopt(cOptDef,varargin{:}); fid       = opt.fid;
  
  
  % Model parameters. rho       = opt.rho;     % Density. miu       = [num2str(opt.mu0*opt.K),'*(sqrt(2)*sqrt(1/2*(4*ux^2+(uy+vx)^2+4*vy^2)+eps))^(',num2str(opt.n-1),')']; % Effective viscosity. % 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&lt;0 )
      fea.grid = quad2tri( fea.grid );   end end n_bdr = max(fea.grid.b(3,:));           % Number of boundaries.
  
  
  % Boundary conditions. dtol      = opt.hmax; i_inflow  = findbdr( fea, ['x&lt;',num2str(dtol)] );     % Inflow boundary number. i_outflow = findbdr( fea, ['x&gt;',num2str(l-dtol)] );   % Outflow boundary number.
  
  
  % Problem definition. fea = addphys(fea,@navierstokes);     % Add Navier-Stokes equations physics mode. fea.phys.ns.eqn.coef{1,end} = { rho }; fea.phys.ns.eqn.coef{2,end} = { miu }; fea.phys.ns.eqn.coef{7,end}{1} = [num2str(opt.dp),'*(1-x/',num2str(l),')'];   % Initial perssure. fea.phys.ns.sfun = { sf_u sf_u sf_p };           % Set shape functions. fea.phys.ns.bdr.sel([i_inflow,i_outflow]) = 4; fea.phys.ns.bdr.coef{4,end}{3,i_inflow}   = l/opt.dp;   % Set inflow pressure.
  
  
  % Parse and solve problem. fea = parsephys(fea);   % Check and parse physics mode. fea = parseprob(fea);   % Check and parse problem struct.
  
  [fea.bdr.n{1}{[i_inflow,i_outflow]}] = deal( '-p*nx' );   % Offset natural Neumann BC. [fea.bdr.d{2}{[i_inflow,i_outflow]}] = deal(0);    % Set v=0 at inlet and outlet.
  
  fea.sol.u = solvestat( fea, 'relchg', false, 'fid', fid );   % Call to stationary solver.
  
  
  % Postprocessing. s_velm = 'sqrt(u^2+v^2)'; s_refsol  = [num2str(opt.n/(opt.n+1)*(1/(opt.mu0*opt.K)*opt.dp)^(1/opt.n)),'*(',num2str((h/2)^((opt.n+1)/opt.n)),'-abs(',num2str(h/2),'-y)^',num2str((opt.n+1)/opt.n),')']; p = [ l/2*ones(1,25); linspace(0,h,25) ]; u = evalexpr( s_velm, p, fea ); u_ref = evalexpr( s_refsol, p, fea ); s_velm = 'sqrt(u^2+v^2)'; if ( opt.iplot&gt;0 )   figure   subplot(3,1,1)   postplot(fea,'surfexpr',s_velm,'evaltype','exact')   title('Velocity field')
  
    subplot(3,1,2)   postplot(fea,'surfexpr','miu_ns')   title('Effective viscosity')
  
    subplot(3,1,3)   plot( u, p(2,:) )   hold on   plot( u_ref, p(2,:), 'k.' )   legend( 'Computed solution', 'Reference solution','Location','West')   title( 'Velocity profile at x=l/2' )   ylabel( 'y' ) end
  
  
  % Error checking. err = sqrt(sum((u-u_ref).^2)/sum(u_ref.^2)); if( ~isempty(fid) )   fprintf(fid,'\nL2 Error: %f\n',err)   fprintf(fid,'\n\n') end
  
  
  out.err  = err; out.pass = err&lt;0.05; if ( nargout==0 )   clear fea out end
  </code></pre>