← Back to branch summary

~fluidity-core/fluidity/exorcised

« back to all changes in this revision

Viewing changes to tools/LocalMats.F90

Latest trunk changes, and modifications to the backward_facing_step_2d example, which are experimental and work is ongoing. It's better than the current trunk example in that at least the k_epsilon code is correct.

Show diffs side-by-side

added added

removed removed

Lines of Context:
tools/LocalMats.F90
1
 
  program LocalMats
2
 
    ! A small program for constructing local matrices
3
 
    ! Useful for doing linear analyses on mixed elements etc.
4
 
    use fields
5
 
    use FEtools
6
 
    use DGtools
7
 
    use elements
8
 
    use transform_elements
9
 
    use signal_vars
10
 
    !use global_parameters, only : current_debug_level, PYTHON_FUNC_LEN
11
 
    !use fldebug
12
 
 
13
 
    implicit none
14
 
    integer :: degree, quad_degree=6
15
 
    type(quadrature_type), target :: quad,f_quad
16
 
    type(element_type), target :: X_shape, u_shape, h_shape
17
 
    real, allocatable, dimension(:,:) :: X_ele
18
 
    character(len=100) :: tmpbuf
19
 
    ! Arguments for handling the command line
20
 
    character(len=256) :: filename
21
 
    integer :: status, ele
22
 
    character(len=100) :: fmt,buffer,buf
23
 
    integer :: u_degree=1, h_degree=2, dim=1, simplex_type=1
24
 
 
25
 
    !simplex types are:
26
 
    !1. Equilateral, unit. 2. Right-angled, unit.
27
 
 
28
 
    namelist/LocalMats_data/u_degree,h_degree,dim,simplex_type, quad_degree
29
 
    logical :: file_exists
30
 
    integer :: unit, io1, stat
31
 
 
32
 
    !current_debug_level = 2
33
 
 
34
 
    print *,  'program LocalMats'
35
 
 
36
 
    call get_command_argument(1, value=filename, status=status)
37
 
    select case(status)
38
 
    case(1:)
39
 
       print *,  'no file specified?'
40
 
       stop
41
 
    case(:-1)
42
 
       write(0,*) "Warning: truncating filename"
43
 
    end select
44
 
    filename=trim(filename)
45
 
 
46
 
    print *,  'Reading parameters'
47
 
 
48
 
    unit = free_unit()
49
 
    open(unit=unit, file=trim(trim(filename)//".dat"), status='old', &
50
 
         iostat=io1)
51
 
 
52
 
    if(io1.ne.0) then
53
 
       print *, 'Looked for ', trim(trim(filename)//".dat")
54
 
       print *, 'Could not read from .dat file'
55
 
       stop
56
 
    end if
57
 
 
58
 
    read(unit, LocalMats_data)
59
 
    close(unit) 
60
 
 
61
 
    print *,  'Making local matrices for an element with'
62
 
    print *,  'Dimension == ', dim
63
 
    print *,  'U degree = ', u_degree
64
 
    print *,  'H degree = ', h_degree
65
 
    print *,  'simplex type = ',simplex_type 
66
 
    print *,  'Quad degree =', quad_degree
67
 
 
68
 
    print *,  'Getting quadrature'
69
 
 
70
 
    quad=make_quadrature(loc=dim+1, dimension=dim, degree=quad_degree)
71
 
 
72
 
    print *,  'Getting shape functions'
73
 
 
74
 
    ! Shape functions for positions (linear)
75
 
    X_shape=make_element_shape(loc=dim+1, dimension=dim, &
76
 
         degree=1, quad=quad)
77
 
 
78
 
    ! Shape functions for velocity and height
79
 
    u_shape=make_element_shape(loc=dim+1, &
80
 
         dimension=dim, degree=u_degree, quad=quad)
81
 
    h_shape=make_element_shape(loc=dim+1, &
82
 
         dimension=dim, degree=h_degree, quad=quad)
83
 
 
84
 
    allocate(X_ele(u_shape%dim,x_shape%loc))
85
 
 
86
 
    call get_X_ele(X_ele,simplex_type)
87
 
 
88
 
    print *,  'Coordinates'
89
 
    print *,  X_ele
90
 
 
91
 
    call assemble_local_matrices(u_shape,h_shape,X_shape,X_ele)
92
 
 
93
 
  contains
94
 
 
95
 
    subroutine assemble_local_matrices( &
96
 
         u_shape,h_shape,X_shape,X_ele)
97
 
      type(element_type), intent(in) :: u_shape,h_shape,X_shape
98
 
      real, dimension(u_shape%dim,x_shape%loc), intent(in) :: X_ele
99
 
 
100
 
      !! Local Stuff
101
 
      ! Coordinate transform * quadrature weights.
102
 
      real, dimension(u_shape%ngi) :: detwei
103
 
      ! Derivatives of shape function:
104
 
      real, dimension(h_shape%loc, h_shape%ngi, h_shape%dim) :: h_dshape
105
 
      ! gradient matrix
106
 
      real, dimension(u_shape%dim,h_shape%loc,u_shape%loc) :: grad_mat
107
 
      real, dimension(u_shape%loc,u_shape%loc) :: u_mass_mat
108
 
      real, dimension(h_shape%loc,h_shape%loc) :: lap_mat,h_mass_mat
109
 
      integer :: idim, iloc
110
 
 
111
 
      ! Transform derivatives and weights into physical space.
112
 
      call transform_to_physical(X_ele, X_shape, m=h_shape, &
113
 
           dm_t=h_dshape, detwei=detwei)
114
 
 
115
 
      grad_mat = -dshape_shape(h_dshape,u_shape,detwei)
116
 
      lap_mat = dshape_dot_dshape(h_dshape,h_dshape,detwei)
117
 
      u_mass_mat = shape_shape(u_shape,u_shape,detwei)
118
 
      h_mass_mat = shape_shape(h_shape,h_shape,detwei)
119
 
 
120
 
      print *,  'U Mass Matrix'
121
 
 
122
 
      do iloc = 1, u_shape%loc
123
 
         print *,  u_mass_mat(iloc,:)
124
 
      end do
125
 
 
126
 
      print *,  'H mass matrix'
127
 
      do iloc = 1, h_shape%loc
128
 
         print *,  h_mass_mat(iloc,:)
129
 
      end do
130
 
 
131
 
      print *,  'Div matrix'
132
 
      do idim = 1, u_shape%dim
133
 
         do iloc = 1, h_shape%loc
134
 
            print *,  grad_mat(idim,iloc,:)
135
 
         end do
136
 
         print *,  ' '
137
 
      end do
138
 
 
139
 
      print *,  'Laplacian matrix'
140
 
      do iloc = 1, h_shape%loc
141
 
         print *, lap_mat(iloc,:)
142
 
      end do
143
 
 
144
 
    end subroutine assemble_local_matrices
145
 
 
146
 
    subroutine get_X_ele(X_ele,simplex_type)
147
 
      real, dimension(:,:), intent(inout) :: X_ele
148
 
      integer, intent(in) :: simplex_type
149
 
 
150
 
      !
151
 
      select case(size(X_ele,1))
152
 
      case (1)
153
 
         X_ele(1,:) = (/ 0.0, 1.0 /)
154
 
      case (2)
155
 
         select case(simplex_type)
156
 
         case (1)
157
 
            X_ele(1,:) = (/ 0.0, 1.0, 0.5 /)
158
 
            X_ele(2,:) = (/ 0.0, 0.0, 0.8660254037844386 /)
159
 
         case (2)
160
 
            X_ele(1,:) = (/ 0.0, 1.0, 1.0 /)
161
 
            X_ele(2,:) = (/ 0.0, 0.0, 1.0 /)
162
 
         case default
163
 
            print *, 'no support for simplex type, try coding it?'
164
 
            stop
165
 
         end select
166
 
      case default
167
 
         print *, 'no support for dimension, try coding it?'
168
 
         stop
169
 
      end select
170
 
 
171
 
    end subroutine get_X_ele
172
 
end program LOCALMATS