merge with beta, QuadRuleDegree set back to 4 for linear tri

BuildMuaWorkers.m

@@ -2,6 +2,7 @@
 
 
 
+
 function MUAworkers=BuildMuaWorkers(CtrlVar,MUA,MUAworkers)
 
 

FE_inner_product.m

@@ -32,21 +32,32 @@
 
 [nM,mM]=size(M);
 
+theta=nan;
+
 if nM==nf
-    
+
     fg=f'*M*g;   % = (M*f)' * g
-    fNorm=sqrt(f'*M*f) ; 
-    gNorm=sqrt(g'*M*g) ; 
-    theta=acos(fg/(fNorm*gNorm)) ; 
-
+
+    if nargout>1
+        fNorm=sqrt(f'*M*f) ;
+        gNorm=sqrt(g'*M*g) ;
+        theta=acos(fg/(fNorm*gNorm)) ;
+    else
+        theta=nan;
+    end
+
 elseif nf==2*nM
-    
+
     fg=f(1:nM)'*M*g(1:nM)+f(1+nM:nf)'*M*g(1+nM:ng);
-
+
+    if nargout>1
     fNorm=sqrt(f(1:nM)'*M*f(1:nM)+f(1+nM:nf)'*M*f(1+nM:nf)) ;
     gNorm=sqrt(g(1:nM)'*M*g(1:nM)+g(1+nM:ng)'*M*g(1+nM:ng)) ;
-    theta=acos(fg/(fNorm*gNorm)) ; 
-
+    theta=acos(fg/(fNorm*gNorm)) ;
+    else
+        theta=nan;
+    end
+
 end
 
 

QuadratureRuleDegree.m

@@ -12,8 +12,8 @@
 
         case 3
 
-            Degree=4;
-            Degree=25;
+            Degree=4;     % 6 integration points
+            % Degree=25;  % this would give 126 integration points
 
 
         case 6

UaUtilities/InsideOutside.m

@@ -6,6 +6,21 @@
 % Just a simple wrapper around inpoly2 to take care of the possibility of
 % several boundaries within the array 'boundary' seperated by NaNs
 %
+%   xy            : n x 2 array
+%   Boundary      : m x 2 Array
+%
+% Example
+% 
+%
+%   xy=[ 0.5 0.5 ; 0.5 1.5] ;
+%   boundary=[ 0 0 ; 1 0 ; 1 1 ; 0 1 ; 0 0 ; nan nan ; 0 2 ; 1 2 ; 1 3 ; 0 3 ; 0 2] ;
+%   [isInside,isOnBounday]=InsideOutside(xy,boundary) ; 
+%
+%   figure ; 
+%   plot(boundary(:,1),boundary(:,2),LineWidth=2) ; hold on ; 
+%   plot(xy(:,1),xy(:,2),"*r")
+%   plot(xy(isInside,1),xy(isInside,2),"ok") ; axis padded
+%
 %%
 
 Kisnan=find(isnan(boundary(:,1)));

UaUtilities/ModifyColormap.m

@@ -24,7 +24,8 @@
 %
 % Use with cmocean, setting colormap to gray around -100 to 100
 %
-%   CM=cmocean('balanced',25,'pivot',0) ; colormap(CM); ModifyColormap(100,nan,ChangeColormap=false) ;
+% CM=cmocean('balanced',25,'pivot',0) ; colormap(CM);
+% ModifyColormap(100,nan,ChangeColormap=false) ;
 %
 %%
 

UaUtilities/PlotCalvingFronts.m

@@ -1,3 +1,7 @@
+
+
+
+
 function [xc,yc]=PlotCalvingFronts(CtrlVar,MUA,F,varargin)
 
 %

UpdateMUA.m

@@ -18,7 +18,7 @@
 % first make sure that the element is of the right type
 [MUA.coordinates,MUA.connectivity]=ChangeElementType(MUA.coordinates,MUA.connectivity,CtrlVar.TriNodes);
 
-% checking if MUA as the QuadratureRuleDegree field
+% checking if MUA has the QuadratureRuleDegree field
 if CtrlVar.QuadRules2021  && ~isfield(MUA,"QuadratureRuleDegree")
     QuadratureFieldMissing=true;
 else