1 | SUBROUTINE addcell(id, eltype, nb_faces, faceslist)
2 |
3 | ! ================================================================!
4 | ! !
5 | ! addcell.F : Adds a new cell to the 'cell_list' and sets the !
6 | ! cell nodal values and calculates the nodal normals !
7 | ! !
8 | ! in : The new cell's 'id' number !
9 | ! The type of cell to construct !
10 | ! The number of faces that compose the cell !
11 | ! The faces' id list !
12 | ! nota : Be careful when using this subroutine. the pointer !
13 | ! 'current_cell' must be pointing at a previously !
14 | ! allocated memory space of type 'cell'. Take a look !
15 | ! at the examples.
16 | ! out : the 'current_cell' and 'last_cell' pointers will !
17 | ! point at the cell !
18 | ! !
19 | ! author : J. AMAYA (avril 2007) !
20 | ! !
21 | ! ================================================================!
22 |
23 | USE datas
24 | IMPLICIT NONE
25 |
26 | INCLUDE 'dom_constants.h'
27 |
28 | ! IN
29 | DOM_INT :: parent, nb_faces, id, eltype
30 | DOM_INT, DIMENSION(nb_faces) :: faceslist
31 |
32 | ! LOCAL
33 | type(cell), pointer :: new_cell
34 | integer :: ierr, nbnodes, ndim
35 | integer :: i, j, k, m, thisnode_id
36 | DOM_REAL :: vol
37 | DOM_REAL, allocatable, dimension(:,:) :: nodescoor
38 | DOM_REAL, allocatable, dimension(:,:) :: normlist
39 | type(face), pointer :: face_pt
40 | logical :: found
41 |
42 | ! ------------------------!
43 | ! Initialise the new cell !
44 | ! ------------------------!
45 |
46 | new_cell => current_cell
47 | new_cell%i_nbfaces = nb_faces
48 | new_cell%cell_id = id
49 | ALLOCATE(new_cell%cell_face(nb_faces))
50 |
51 | ! ----------------------!
52 | ! Add faces to the cell !
53 | ! ----------------------!
54 |
55 | ! print*, " >> cell: ", id
56 | ! print*, " faces: ", faceslist
57 |
58 | DO i=1, nb_faces
59 | CALL findface(faceslist(i), ierr)
60 | ! print*, " face ",i," to add: ", faceslist(i)
61 | IF (ierr.eq.0) THEN
62 | WRITE(*,*) "Error in addcell:face not found:",faceslist(i)
63 | ELSE
64 | new_cell%cell_face(i)%face_ptr => current_face
65 | ENDIF
66 | ENDDO
67 |
68 | ! ---------------!
69 | ! Set cell nodes !
70 | ! ---------------!
71 | new_cell%celltype = eltype
72 |
73 | SELECTCASE(eltype)
74 | CASE(EL_TRI)
75 | ndim = 2
76 | nbnodes = 3
77 | CASE(EL_QUAD)
78 | ndim = 2
79 | nbnodes = 4
80 | CASE(EL_TETRA)
81 | ndim = 3
82 | nbnodes = 4
83 | CASE(EL_HEXA)
84 | ndim = 3
85 | nbnodes = 8
86 | CASE DEFAULT
87 | WRITE(*,*) "Error: element type not recognised ", eltype
88 | STOP
89 | ENDSELECT
90 |
91 | ! print*, "Allocating cellnodes: ", nbnodes
92 | ALLOCATE(new_cell%cellnodes(nbnodes))
93 | ALLOCATE(new_cell%nodenormals(4,nbnodes))
94 |
95 | new_cell%cellnodes = 0
96 |
97 | k = 1
98 | i = 1
99 |
100 | DO WHILE ((i.le.nb_faces).and.(k.le.nbnodes))
101 |
102 | face_pt => new_cell%cell_face(i)%face_ptr
103 | j = 1
104 |
105 | DO WHILE ((j.le.face_pt%i_nbnodes).and.(k.le.nbnodes))
106 |
107 | parent = 1
108 | IF (face_pt%parent_cell(1).ne.id) parent = 2
109 |
110 | ! print*, " >> Evaluating: face ", i,", node ",j
111 | ! print*, " parent: " ,parent, " , k: ", k
112 | thisnode_id = face_pt%face_point(j,parent)
113 | found = .false.
114 | m = k - 1
115 | ! print*, " m: ",m
116 | DO WHILE ((m.ge.1).and.(.not.found))
117 | IF (thisnode_id.eq.new_cell%cellnodes(m)) found = .true.
118 | m = m - 1
119 | ENDDO
120 | IF (.not.found) THEN
121 | ! print*, " adding to the list bcs not found"
122 | new_cell%cellnodes(k) = thisnode_id
123 | k = k + 1
124 | ENDIF
125 | ! print*, " done"
126 | ! print*, " "
127 | j = j + 1
128 |
129 | ENDDO
130 |
131 | i = i + 1
132 |
133 | ENDDO
134 |
135 | ! ----------------------------!
136 | ! Set cell's node information !
137 | ! ----------------------------!
138 |
139 | ALLOCATE(nodescoor(ndim,nbnodes))
140 | ALLOCATE(normlist(ndim,nbnodes))
141 |
142 | ! -------------------------!
143 | ! Calculate normal at node !
144 | ! -------------------------!
145 |
146 | DO i=1, nbnodes
147 | thisnode_id = new_cell%cellnodes(i)
148 | ! print*, " node ",i,": ", thisnode_id
149 | nodescoor(:,i) = node_list(:,thisnode_id)
150 |
151 | new_cell%nodenormals(:,i) = 0.0d0
152 | DO j=1, nb_faces
153 |
154 | face_pt => new_cell%cell_face(j)%face_ptr
155 |
156 | parent = 1
157 | IF (face_pt%parent_cell(1).ne.id) parent = 2
158 | ! print*, " faces parents:", face_pt%parent_cell(1), &
159 | ! & face_pt%parent_cell(2)
160 | ! print*, " parent: ",parent
161 |
162 | ! --------------------------------------!
163 | ! Detect if this face contains the node !
164 | ! --------------------------------------!
165 | ierr = 0
166 | k=1
167 | DO WHILE (k.le.face_pt%i_nbnodes)
168 | IF (face_pt%face_point(k,parent).eq.thisnode_id) THEN
169 | ierr=thisnode_id
170 | k = face_pt%i_nbnodes
171 | ENDIF
172 | k = k+1
173 | ENDDO
174 |
175 | ! ---------------------------------------------------------!
176 | ! Add current face's area weighted normal to node's normal !
177 | ! ---------------------------------------------------------!
178 | IF (ierr.ne.0) THEN
179 |
180 | new_cell%nodenormals(1,i) = new_cell%nodenormals(1,i) + &
181 | & face_pt%face_normal(1,parent)*face_pt%area
182 | new_cell%nodenormals(2,i) = new_cell%nodenormals(2,i) + &
183 | & face_pt%face_normal(2,parent)*face_pt%area
184 | new_cell%nodenormals(3,i) = new_cell%nodenormals(3,i) + &
185 | & face_pt%face_normal(3,parent)*face_pt%area
186 |
187 | ENDIF
188 |
189 | ENDDO
190 |
191 | ! ------------------------------------!
192 | ! Calculate node's normal vector norm !
193 | ! ------------------------------------!
194 | new_cell%nodenormals(4,i) = SQRT(new_cell%nodenormals(1,i)* &
195 | & new_cell%nodenormals(1,i)+ &
196 | & new_cell%nodenormals(2,i)* &
197 | & new_cell%nodenormals(2,i)+ &
198 | & new_cell%nodenormals(3,i)* &
199 | & new_cell%nodenormals(3,i))
200 |
201 | ! print*, " nx: ", new_cell%nodenormals(1,i)
202 | ! print*, " ny: ", new_cell%nodenormals(2,i)
203 | ! print*, " nz: ", new_cell%nodenormals(3,i)
204 |
205 | normlist(1:3,i) = new_cell%nodenormals(1:3,i)
206 | ENDDO
207 |
208 | ! ------------------------!
209 | ! Calculate cell's volume !
210 | ! ------------------------!
211 |
212 | CALL calculatevol(ndim,eltype,nbnodes,nodescoor,normlist,vol)
213 | new_cell%volume = vol
214 | print*, " >> Cell ", id, " - volume: ", vol
215 |
216 | DEALLOCATE(nodescoor)
217 | DEALLOCATE(normlist)
218 |
219 | ! -----------------------------!
220 | ! Put the new cell in the list !
221 | ! (This can provoque memory errors) !
222 | ! -----------------------------!
223 |
224 | ! IF (.not.ASSOCIATED(cell_list)) THEN
225 | ! cell_list => new_cell
226 | ! last_cell => new_cell
227 | ! ELSE
228 | ! last_cell%next_cell => new_cell
229 | ! last_cell => new_cell
230 | ! ENDIF
231 |
232 | ! -----------------------------!
233 | ! Put the new cell in the list !
234 | ! -----------------------------!
235 |
236 | IF (i_ncells.eq.0) THEN
237 | cell_list => new_cell
238 | last_cell => new_cell
239 | ELSE
240 | last_cell%next_cell => new_cell
241 | last_cell => new_cell
242 | ENDIF
243 |
244 | i_ncells = i_ncells + 1
245 |
246 | END SUBROUTINE addcell