1!/*****************************************************************************/ 2! * 3! * Elmer/Ice, a glaciological add-on to Elmer 4! * http://elmerice.elmerfem.org 5! * 6! * 7! * This program is free software; you can redistribute it and/or 8! * modify it under the terms of the GNU General Public License 9! * as published by the Free Software Foundation; either version 2 10! * of the License, or (at your option) any later version. 11! * 12! * This program is distributed in the hope that it will be useful, 13! * but WITHOUT ANY WARRANTY; without even the implied warranty of 14! * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15! * GNU General Public License for more details. 16! * 17! * You should have received a copy of the GNU General Public License 18! * along with this program (in file fem/GPL-2); if not, write to the 19! * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 20! * Boston, MA 02110-1301, USA. 21! * 22! *****************************************************************************/ 23! ****************************************************************************** 24! * 25! * Author: F. Gillet-Chaulet (IGE) 26! * Email: fabien.gillet-chaulet@univ-grenoble-alpes.fr 27! * Web: http://elmerice.elmerfem.org 28! * 29! * Original Date: 01/02/2018 30! * 31! * Find conneceted areas in a mesh 32! ***************************************************************************** 33!!! 34 SUBROUTINE GetConnectedAreas( Model,Solver,dt,TransientSimulation) 35 USE DefUtils 36 IMPLICIT NONE 37 38 TYPE(Model_t) :: Model 39 TYPE(Solver_t), TARGET :: Solver 40 REAL(KIND=dp) :: dt 41 LOGICAL :: TransientSimulation 42 43 TYPE Queue_t 44 INTEGER :: maxsize 45 INTEGER :: top 46 INTEGER,ALLOCATABLE :: items(:) 47 END TYPE Queue_t 48 49 TYPE(Variable_t), POINTER :: Var,Area 50 TYPE(Mesh_t), POINTER :: Mesh 51 TYPE(Queue_t) :: Queue 52 TYPE(Element_t),POINTER :: Element,Parent 53 TYPE(Element_t),POINTER :: Faces(:),Face 54 TYPE(GaussIntegrationPoints_t) :: IP 55 TYPE(Nodes_t),SAVE :: Nodes 56 TYPE(ValueList_t), POINTER :: SolverParams 57 REAL(KIND=dp), ALLOCATABLE :: Basis(:), dBasisdx(:,:), & 58 ddBasisddx(:,:,:) 59 REAL(KIND=dp) :: s, detJ 60 REAL(KIND=dp) :: smallest,largest 61 REAL(KIND=dp), ALLOCATABLE :: RegionArea(:) 62 INTEGER :: EIndex 63 INTEGER :: label 64 INTEGER :: region 65 INTEGER,ALLOCATABLE :: ElementLabel(:) 66 INTEGER :: t,i,n,p 67 INTEGER :: nfaces 68 INTEGER :: NOFActive 69 70 LOGICAL :: stat,Found 71 LOGICAL :: SAVE_REGIONS 72 73 CHARACTER(LEN=MAX_NAME_LEN) :: SolverName='GetConnectedAreas' 74 75 Mesh => Solver % Mesh 76 SolverParams => GetSolverParams() 77 78 n = MAX(Mesh % MaxElementNodes,Mesh % MaxElementDOFs) 79 ALLOCATE(ElementLabel(Solver%Mesh%NumberOfBulkElements)) 80 ALLOCATE( Basis(n), dBasisdx(n,3), ddBasisddx(n,3,3)) 81 ALLOCATE(Nodes%x(n),Nodes%y(n),Nodes%z(n)) 82 83 84 CALL FindMeshEdges(Mesh,.FALSE.) 85 SELECT CASE(Mesh % MeshDim) 86 CASE(2) 87 Faces => Mesh % Edges 88 CASE(3) 89 Faces => Mesh % Faces 90 END SELECT 91 92 CALL QueueInit(Queue,Mesh%NumberOfBulkElements) 93 94 NOFActive=GetNOFActive() 95 96 ElementLabel=-1 97 label=1 98 99 DO t=1,NOFActive 100 Element => GetActiveElement(t) 101 EIndex= Element % ElementIndex 102 103 IF (CheckPassiveElement(Element)) CYCLE 104 IF (ElementLabel(EIndex).GT.0) CYCLE 105 106 ElementLabel(EIndex) = label 107 CALL QueuePush(Queue,EIndex) 108 109 DO WHILE (Queue%top.GT.0) 110 CALL QueuePop(Queue,EIndex) 111 Element => Mesh % Elements(EIndex) 112 113 SELECT CASE(Mesh % MeshDim) 114 CASE(2) 115 nfaces=Element % TYPE % NumberOfEdges 116 CASE(3) 117 nfaces=Element % TYPE % NumberOfFaces 118 END SELECT 119 120 DO i=1,nfaces 121 SELECT CASE(Mesh % MeshDim) 122 CASE(2) 123 Face => Mesh % Edges (Element % EdgeIndexes(i)) 124 CASE(3) 125 Face => Mesh % Faces (Element % FaceIndexes(i)) 126 END SELECT 127 128 IF (.NOT.ASSOCIATED(Face)) & 129 CALL FATAL(SolverName,'Face not found') 130 131 Parent => Face % BoundaryInfo % Left 132 IF (ASSOCIATED(Parent)) THEN 133 IF ((.NOT.CheckPassiveElement(Parent))& 134 .AND.(ElementLabel(Parent%ElementIndex).LE.0)) THEN 135 ElementLabel(Parent%ElementIndex)=label 136 CALL QueuePush(Queue,Parent%ElementIndex) 137 END IF 138 END IF 139 Parent => Face % BoundaryInfo % Right 140 IF (ASSOCIATED(Parent)) THEN 141 IF ((.NOT.CheckPassiveElement(Parent))& 142 .AND.(ElementLabel(Parent%ElementIndex).LE.0)) THEN 143 ElementLabel(Parent%ElementIndex)=label 144 CALL QueuePush(Queue,Parent%ElementIndex) 145 END IF 146 END IF 147 148 END DO 149 150 END DO 151 152 label = label + 1 153 END DO 154 155 label = label - 1 156 WRITE(Message,'(A,i0,A)') 'There is ',label,' unconnected areas' 157 CALL INFO(SolverName,TRIM(Message),level=3) 158 159 Var => VariableGet( Mesh % Variables,'AreaNumber',UnfoundFatal=.True.) 160 Area => VariableGet( Mesh % Variables,'RegionArea',UnfoundFatal=.True.) 161 162 ALLOCATE(RegionArea(label)) 163 164 RegionArea=0._dp 165 Var % Values = -1 166 Area % Values = 0. 167 168 DO t=1,NOFActive 169 Element => GetActiveElement(t) 170 EIndex= Element % ElementIndex 171 n = GetElementNOFNodes(Element) 172 173 region=ElementLabel(EIndex) 174 IF (region.LE.0) CYCLE 175 176 Var % Values(Var % Perm(Element % NodeIndexes(1:n)))= & 177 region 178 179 Nodes % x(1:n) = Mesh % Nodes % x(Element % NodeIndexes(1:n)) 180 Nodes % y(1:n) = Mesh % Nodes % y(Element % NodeIndexes(1:n)) 181 Nodes % z(1:n) = Mesh % Nodes % z(Element % NodeIndexes(1:n)) 182 IP = GaussPoints( Element ) 183 DO p = 1, IP % n 184 stat = ElementInfo( Element, Nodes, IP % U(p), IP % V(p), & 185 IP % W(p), detJ, Basis, dBasisdx, ddBasisddx, .FALSE.) 186 s = detJ * IP % S(p) 187 188 RegionArea(region) = RegionArea(region) + s 189 END DO 190 191 192 END DO 193 194 !$omp parallel do & 195 !$omp& private(region) 196 DO t=1,Mesh % NumberOfNodes 197 region=NINT(Var % Values(Var % Perm(t))) 198 IF (region.GT.0) THEN 199 Area % Values(Area % Perm(t))=& 200 MAX(RegionArea(region),Area % Values(Area % Perm(t))) 201 END IF 202 END DO 203 !$omp end parallel do 204 205 largest=MAXVAL(RegionArea) 206 smallest=MINVAL(RegionArea) 207 WRITE(Message,'(A,e15.7)') 'largest area :',largest 208 CALL INFO(SolverName,TRIM(Message),level=3) 209 WRITE(Message,'(A,e15.7)') 'smallest area :',smallest 210 CALL INFO(SolverName,TRIM(Message),level=3) 211 212 SAVE_REGIONS=ListGetLogical(SolverParams,'Save regions labels',Found) 213 IF (.NOT.Found) SAVE_REGIONS=.FALSE. 214 IF (SAVE_REGIONS) THEN 215 Open(12,file='regions.txt') 216 DO i=1,label 217 write(12,*) i,RegionArea(i) 218 END DO 219 close(12) 220 END IF 221 222 DEALLOCATE(ElementLabel) 223 DEALLOCATE(RegionArea) 224 DEALLOCATE(Basis,dBasisdx,ddBasisddx) 225 DEALLOCATE(Queue%items) 226 DEALLOCATE(Nodes%x,Nodes%y,Nodes%z) 227 228 229 CONTAINS 230 SUBROUTINE QueueInit(Queue,n) 231 IMPLICIT NONE 232 TYPE(Queue_t) :: Queue 233 INTEGER :: n 234 ALLOCATE(Queue%items(n)) 235 Queue%top=0 236 Queue%maxsize=n 237 END SUBROUTINE QueueInit 238 239 SUBROUTINE QueuePush(Queue,x) 240 IMPLICIT NONE 241 TYPE(Queue_t) :: Queue 242 INTEGER :: x 243 IF (Queue%top.EQ.Queue%maxsize) & 244 CALL FATAL(SolverName,'Too many elements in the queue?') 245 246 Queue%top=Queue%top+1 247 Queue%items(Queue%top)=x 248 249 END SUBROUTINE QueuePush 250 251 SUBROUTINE QueuePop(Queue,x) 252 IMPLICIT NONE 253 TYPE(Queue_t) :: Queue 254 INTEGER :: x 255 IF (Queue%top.EQ.0) & 256 CALL FATAL(SolverName,'No more element in the queue') 257 258 x=Queue%items(Queue%top) 259 Queue%top=Queue%top-1 260 END SUBROUTINE QueuePop 261 END 262 263 264 265