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