automated repeated mesh refinement

AdaptMesh.m

@@ -6,9 +6,6 @@
 narginchk(10,10)
 nargoutchk(7,7)
 
-persistent FigMesh VideoMesh
-
-
 MUAnew=MUAold;
 Fnew=Fold;
 BCsNew=BCsOld;
@@ -20,12 +17,6 @@
     fprintf('Before remeshing: '); PrintInfoAboutElementsSizes(CtrlVar,MUAold)
 end
 
-if CtrlVar.doplots && CtrlVar.doAdaptMeshPlots && CtrlVar.InfoLevelAdaptiveMeshing>=2
-    figure
-    PlotFEmesh(MUAold.coordinates,MUAold.connectivity,CtrlVar);
-    title(sprintf('Initial mesh on input to AdaptMesh : #Ele=%-i, #Nodes=%-i, #nod=%-i',MUAold.Nele,MUAold.Nnodes,MUAold.nod))
-end
-
 
 isMeshAdvanceRetreat = CtrlVar.FEmeshAdvanceRetreat && ( ReminderFraction(CtrlVar.time,CtrlVar.FEmeshAdvanceRetreatDT)<1e-5 || CtrlVar.FEmeshAdvanceRetreatDT==0);
 
@@ -35,6 +26,11 @@
     && ~isMeshAdvanceRetreat;
 
 
+JJ=0 ;
+nNewElements=inf;
+nNewNodes=inf;
+nNoChange=0;
+
 
 if isMeshAdvanceRetreat
 
@@ -47,105 +43,112 @@
 
 elseif isMeshAdapt
 
-    for JJ=1:CtrlVar.AdaptMeshIterations
+    while true
+
+
+        JJ=JJ+1;
+
+        if JJ>=CtrlVar.AdaptMeshIterations
+            fprintf('Breaking out of adapt mesh iteration because number of iterations (%i%) exceeds ''CtrlVar.AdaptMeshIterations'' (%i%)\n',...
+                JJ,CtrlVar.AdaptMeshIterations)
+            break
+        end
+
+        if abs(nNewElements)<CtrlVar.AdaptMeshUntilChangeInNumberOfElementsLessThan
+
+            nNoChange=nNoChange+1;
+            %if ~contains(RunInfo.MeshAdapt,'Bisection Coarsening')
+            if nNoChange>2
+
+                fprintf('Breaking out of adapt mesh iteration because change in the number of elements (%i) less than ''CtrlVar.AdaptMeshUntilChangeInNumberOfElementsLessThan'' (%i)\n',...
+                    abs(nNewElements),CtrlVar.AdaptMeshUntilChangeInNumberOfElementsLessThan)
+                break
+
+            end
+        else
+            nNoChange=0;
+        end
+
+
 
         MUAold=MUAnew;
         Fold=Fnew;
         BCsOld=BCsNew;
         GFold=GFnew;
 
         if CtrlVar.InfoLevelAdaptiveMeshing>=1
-            fprintf(CtrlVar.fidlog,' =====  Remeshing at start of run step %-i. Iteration #%-i out of %-i \n ',CtrlVar.CurrentRunStepNumber,JJ,CtrlVar.AdaptMeshIterations);
+            fprintf(CtrlVar.fidlog,' =====  Remeshing at start of run step %-i. Remeshing iteration #%-i \n ',CtrlVar.CurrentRunStepNumber,JJ);
         end
-
+        
         %  Determine new desired element sizes and identify elements for refinement
-        %  or coarsening. 
+        %  or coarsening.
         [UserVar,RunInfo,Fold,xNod,yNod,EleSizeDesired,ElementsToBeRefined,ElementsToBeCoarsened]=...
             NewDesiredEleSizesAndElementsToRefineOrCoarsen2(UserVar,RunInfo,CtrlVar,MUAold,BCsOld,Fold,lold,GFold,RuvOld,Lubvb);
 
         %  Remesh: either global-remeshing or local-mesh refinement.
         [UserVar,RunInfo,CtrlVar,MUAnew]=...
             Remeshing(UserVar,RunInfo,CtrlVar,MUAold,BCsOld,Fold,lold,GFold,...
             xNod,yNod,EleSizeDesired,ElementsToBeRefined,ElementsToBeCoarsened);
-
+        
 
         if MUAnew.Nele==0
             fprintf('No elements left in mesh! \n ')
             return
         end
 
 
+        nNewElements=MUAnew.Nele-MUAold.Nele;
+        nNewNodes=MUAnew.Nnodes-MUAold.Nnodes;
+
         [UserVar,Fnew,BCsNew,GFnew,lnew]=MapFbetweenMeshes(UserVar,CtrlVar,MUAold,MUAnew,Fold,BCsOld,GFold,lold);
 
+
         %% Plots
-        if (CtrlVar.doplots && CtrlVar.doAdaptMeshPlots && CtrlVar.InfoLevelAdaptiveMeshing>=1) || CtrlVar.CreateMeshAdaptVideo
-
+        if  CtrlVar.doplots && CtrlVar.doAdaptMeshPlots && CtrlVar.InfoLevelAdaptiveMeshing>=10
 
-            if isempty(FigMesh)
-                FigMesh=figure ;
-                %FigMesh.Position=[0 0 figsWidth 3*figHeights];
-                FigMesh.Position=[1 1 1000 1000];% full laptop window
-                axis tight
-                if CtrlVar.CreateMeshAdaptVideo
-                    VideoMesh=VideoWriter('MeshAdapt.avi');
-                    VideoMesh.FrameRate=1;
-                    open(VideoMesh);
-                end
+            FigNameBA='Adapt Mesh: before and after';
+            fig=findobj(0,'name',FigNameBA);
+            if isempty(fig)
+                fig=figure('name',FigNameBA);
+                fig.Position=[100,100,1000,1000] ;
             else
-                FigMesh=figure ;
+                fig=figure(fig);
                 hold off
             end
 
-            if CtrlVar.PlotBCs
-                PlotBoundaryConditions(CtrlVar,MUAnew,BCsNew);
-            else
-
-                PlotFEmesh(MUAnew.coordinates,MUAnew.connectivity,CtrlVar); hold on
-                title(sprintf('AdaptMesh Iteration %i: #Ele=%-i, #Nodes=%-i, #nod=%-i',JJ,MUAnew.Nele,MUAnew.Nnodes,MUAnew.nod))
-            end
-
-            if CtrlVar.GLmeshing
-                hold on ; plot(xGLmesh,yGLmesh,'-rx','LineWidth',2);
-            end
+            subplot(2,1,1)
+            hold off
+            xGL=[] ; yGL=[]; GLgeo=[];
+            PlotMuaMesh(CtrlVar,MUAold,[],CtrlVar.MeshColor);
+            hold on ;  [xGL,yGL]=PlotGroundingLines(CtrlVar,MUAold,GFold,GLgeo,xGL,yGL,'r');
+            title(sprintf('Mesh before remeshing  \t #Ele=%-i, #Nodes=%-i, #nod=%-i',MUAold.Nele,MUAold.Nnodes,MUAold.nod))
+            axis tight
 
+            subplot(2,1,2)
+            hold off
+            xGL=[] ; yGL=[]; GLgeo=[];
+            CtrlVar.PlotGLs=1;
+            PlotMuaMesh(CtrlVar,MUAnew,[],CtrlVar.MeshColor);
+            title(sprintf('Mesh after remeshing  \t #Ele=%-i, #Nodes=%-i, #nod=%-i',MUAnew.Nele,MUAnew.Nnodes,MUAnew.nod))
+            hold on ;  [xGL,yGL]=PlotGroundingLines(CtrlVar,MUAnew,GFnew,GLgeo,xGL,yGL,'r');
             axis tight
-            if CtrlVar.CreateMeshAdaptVideo
-                frame = getframe(gcf);
-                writeVideo(VideoMesh,frame);
-            end
+
+
+
         end
         %%
-
     end
-
 end
 
 
-
 if CtrlVar.InfoLevelAdaptiveMeshing>=1
-    fprintf('After remeshing: ') ; 
+    fprintf('After remeshing: ') ;
     PrintInfoAboutElementsSizes(CtrlVar,MUAnew)
 end
 
 
 if CtrlVar.AdaptMeshAndThenStop
 
-    if CtrlVar.doplots  && CtrlVar.InfoLevelAdaptiveMeshing>=10
-
-
-        xGL=[] ; yGL=[]; GLgeo=[];
-        CtrlVar.PlotGLs=1;
-        figure ; PlotMuaMesh(CtrlVar,MUAnew,[],CtrlVar.MeshColor);
-        title(sprintf('Mesh after remeshing  \t #Ele=%-i, #Nodes=%-i, #nod=%-i',MUAnew.Nele,MUAnew.Nnodes,MUAnew.nod))
-        hold on ;  [xGL,yGL]=PlotGroundingLines(CtrlVar,MUAnew,GFnew,GLgeo,xGL,yGL,'r');
-
-        xGL=[] ; yGL=[]; GLgeo=[];
-        figure ; PlotMuaMesh(CtrlVar,MUAold,[],CtrlVar.MeshColor);
-        hold on ;  [xGL,yGL]=PlotGroundingLines(CtrlVar,MUAold,GFold,GLgeo,xGL,yGL,'r');
-        title(sprintf('Mesh before remeshing  \t #Ele=%-i, #Nodes=%-i, #nod=%-i',MUAold.Nele,MUAold.Nnodes,MUAold.nod))
-        %%
-    end
-
     MUA=MUAnew;
     save(CtrlVar.SaveInitialMeshFileName,'MUA') ;
     fprintf(CtrlVar.fidlog,'New mesh was saved in %s .\n',CtrlVar.SaveAdaptMeshFileName);
@@ -161,7 +164,7 @@
 
 %  Do velocities need to be recalculated?
 %
-%  Always recalculate velocities if: 
+%  Always recalculate velocities if:
 %
 %   CtrlVar.InitialDiagnosticStepAfterRemeshing is true
 %   but also if mesh refinement method was not 'newest vertex bisection'
@@ -181,10 +184,6 @@
     fprintf(' Adapted FE mesh was saved in %s .\n',CtrlVar.SaveAdaptMeshFileName);
 end
 
-if ~isempty(VideoMesh)
-    close(VideoMesh)
-end
-
 
 end
 

Error2EleSize.m

@@ -29,7 +29,7 @@
 % for p=2 : h(e0) =  1/4 (hMax-hMin)  + hMin
 %
 %
-%
+%%
 
 if isempty(hMin)
     hMin=CtrlVar.MeshSizeMin;

LocalMeshRefinement.m

@@ -1,4 +1,4 @@
-function MUAnew=LocalMeshRefinement(CtrlVar,MUAold,ElementsToBeRefined,ElementsToBeCoarsened)
+function [MUAnew,RunInfo]=LocalMeshRefinement(CtrlVar,RunInfo,MUAold,ElementsToBeRefined,ElementsToBeCoarsened)
 
 persistent wasRefine
 
@@ -31,7 +31,8 @@
 % CtrlVar.LocalMeshRefinementMethod='red-green';
 % CtrlVar.LocalMeshRefinementMethod='newest vertex bisection';
 
-narginchk(4,4)
+narginchk(5,5)
+
 
 if ~islogical(ElementsToBeRefined)
 
@@ -58,15 +59,14 @@
 end
 
 
-%figure ; subplot(1,2,1) ; PlotMuaMesh(CtrlVar,UpdateMUA(CtrlVar,MUAold));  title('original')
-
 [MUAold.coordinates,MUAold.connectivity]=ChangeElementType(MUAold.coordinates,MUAold.connectivity,3);
 
 switch CtrlVar.MeshRefinementMethod
 
     case {'explicit:local:red-green','explicit:local'}
 
 
+        RunInfo.MeshAdapt='Red-Green Refinement';
         [MUAold.coordinates,MUAold.connectivity] = refine(MUAold.coordinates,MUAold.connectivity,T);
         [MUAold.coordinates] = GHGsmoothmesh(MUAold.coordinates,MUAold.connectivity,CtrlVar.LocalAdaptMeshSmoothingIterations,[]);
         MUAold.connectivity=FlipElements(MUAold.connectivity);
@@ -96,26 +96,21 @@
 
         Na=size(mesh.coordinates,1);  Ea=size(mesh.elements,1);
 
-        % do refinement/coarsening alternativily
-        if wasRefine
-            isRefine=0;
-        else
-            isRefine=1;
-        end
-
-        % except if either nRefine or nCoarsen happens to be zero
+        % do refinement/coarsening alternativily, unless if there are
+        % significant in number of elements to be refined/coarsened
 
-        if nCoarsen==0
-            isRefine=1;
-        end
 
-        if nRefine==0
-            isRefine=0;
+        if wasRefine  &&  ( nRefine < 2*nCoarsen) || nRefine==0
+                isRefine=0;
+        else
+                isRefine=1;
         end
+
 
         wasRefine=isRefine;
 
         if isRefine
+            RunInfo.MeshAdapt='Bisection Refinement';
             fprintf(' Refining %i elements \n',nRefine)
             if MUAold.nod~=3
                 meshElementsToBeRefined=pointLocation(mesh.TR,[MUAold.xEle(ElementsToBeRefined) MUAold.yEle(ElementsToBeRefined)]);
@@ -124,6 +119,7 @@
             end
             mesh = bisectionRefine2D(mesh,meshElementsToBeRefined);
         else
+            RunInfo.MeshAdapt='Bisection Coarsening';
             fprintf(' Coarsening %i elements \n',nCoarsen)
             %mesh.elements=FlipElements(mesh.elements);
             if MUAold.nod~=3
@@ -159,12 +155,5 @@
 
 end
 
-%mesh1 = bisectionCoarsen(mesh1,'all') ;  %markedElements, varargin)
-%figure(1000) ; PlotFEmesh(mesh1.coordinates,mesh1.elements,CtrlVar) ;  title('mesh coarsen')
-
-%subplot(1,2,2) ; PlotMuaMesh(CtrlVar,UpdateMUA(CtrlVar,MUAnew));  title('new')
-%%
-
-
 
 end

MapFbetweenMeshes.m

@@ -28,8 +28,8 @@
     if CtrlVar.CurrentRunStepNumber==1  && ~CtrlVar.Restart
 
         fprintf('Note: As this is the first run-step in a time-dependent run: \n')
-        fprintf('When mapping quantities from an old to a new mesh, all geometrical variables (s, b, S, and B) of the new mesh \n')
-        fprintf('are therefore obtained through a call to DefineGeometry.m and not through interpolation from the old mesh.\n')
+        fprintf('        When mapping quantities from an old to a new mesh, all geometrical variables (s, b, S, and B) of the new mesh \n')
+        fprintf('        are defined through a call to DefineGeometry.m and not through interpolation from the old mesh.\n')
         [UserVar,Fnew.s,Fnew.b,Fnew.S,Fnew.B,Fnew.alpha]=GetGeometry(UserVar,CtrlVar,MUAnew,CtrlVar.time,'sbSB');
         Fnew.h=Fnew.s-Fnew.b;
 
@@ -44,8 +44,8 @@
 
     % if a diagnostic step then surface (s) and bed (b), and hence the thickness (h), are defined by the user
     fprintf('Note: As this is not a time-dependent run: \n')
-    fprintf('When mapping quantities from an old to a new mesh, all geometrical variables (s, b, S, and B) of the new mesh \n')
-    fprintf('are therefore obtained through a call to DefineGeometry.m and not through interpolation from the old mesh.\n')
+    fprintf('        When mapping quantities from an old to a new mesh, all geometrical variables (s, b, S, and B) of the new mesh \n')
+    fprintf('        are defined through a call to DefineGeometry.m and not through interpolation from the old mesh.\n')
 
     [UserVar,Fnew.s,Fnew.b,Fnew.S,Fnew.B,Fnew.alpha]=GetGeometry(UserVar,CtrlVar,MUAnew,CtrlVar.time,'sbSB');
     Fnew.h=Fnew.s-Fnew.b;