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