FrontISTR  5.2.0
Large-scale structural analysis program with finit element method
set_arrays_DirectSolver.f90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
2 ! Copyright (c) 2019 FrontISTR Commons
3 ! This software is released under the MIT License, see LICENSE.txt
4 !-------------------------------------------------------------------------------
7 
9 
10  use m_fstr
12 
13  implicit none
14 
15  integer (kind=kint) :: numnon0
16  integer (kind=kint), allocatable :: pointers(:)
17  integer (kind=kint), allocatable :: indices(:)
18  real (kind=kreal) , allocatable :: values(:)
19 
21 
22 contains
23 
24 
27  subroutine set_pointersandindices_directsolver(hecMAT,fstrMAT,is_sym)
28 
29  type(hecmwst_matrix) :: hecMAT
30  type (fstrST_matrix_contact_lagrange) :: fstrMAT
31  logical :: is_sym
32 
33  integer (kind=kint) :: np
34  integer (kind=kint) :: ndof
35  integer (kind=kint) :: num_lagrange
36  integer (kind=kint) :: nn
37  integer (kind=kint) :: ierr
38  integer (kind=kint) :: i, j, k, l, countNon0
39 
40  symmetricmatrixstruc = is_sym
41 
42  np = hecmat%NP ; ndof = hecmat%NDOF ; num_lagrange = fstrmat%num_lagrange
43  nn = np*ndof + num_lagrange + 1
44 
45  if( symmetricmatrixstruc )then
46  numnon0 = hecmat%NPU*ndof**2+hecmat%NP*ndof*(ndof+1)/2 &
47  + (fstrmat%numU_lagrange)*ndof+fstrmat%num_lagrange
48  else
49  numnon0 = (hecmat%NPL+hecmat%NPU+hecmat%NP)*ndof**2 &
50  + (fstrmat%numL_lagrange+fstrmat%numU_lagrange)*ndof
51  endif
52 
53  if(allocated(pointers))deallocate(pointers)
54  allocate(pointers(nn), stat=ierr)
55  if( ierr /= 0 ) stop " Allocation error, mkl%pointers "
56  pointers = 0
57 
58  if(allocated(indices))deallocate(indices)
59  allocate(indices(numnon0), stat=ierr)
60  if( ierr /= 0 ) stop " Allocation error, mkl%indices "
61  indices = 0
62 
63  pointers(1) = 1
64  countnon0 = 1
65 
66  do i = 1, np
67  do j = 1, ndof
68  if( .not. symmetricmatrixstruc )then
69  do l = hecmat%indexL(i-1)+1, hecmat%indexL(i)
70  do k = 1, ndof
71  indices(countnon0) = (hecmat%itemL(l)-1)*ndof + k
72  countnon0 = countnon0 + 1
73  enddo
74  enddo
75  do k = 1, j-1
76  indices(countnon0) = (i-1)*ndof + k
77  countnon0 = countnon0 + 1
78  enddo
79  endif
80  do k = j, ndof
81  indices(countnon0) = (i-1)*ndof + k
82  countnon0 = countnon0 + 1
83  enddo
84  do l = hecmat%indexU(i-1)+1, hecmat%indexU(i)
85  do k = 1, ndof
86  indices(countnon0) = (hecmat%itemU(l)-1)*ndof + k
87  countnon0 = countnon0 + 1
88  enddo
89  enddo
90  if( num_lagrange > 0 )then
91  do l = fstrmat%indexU_lagrange(i-1)+1, fstrmat%indexU_lagrange(i)
92  indices(countnon0) = np*ndof + fstrmat%itemU_lagrange(l)
93  countnon0 = countnon0 + 1
94  enddo
95  endif
96  pointers((i-1)*ndof+j+1) = countnon0
97  enddo
98  enddo
99 
100  if( num_lagrange > 0 )then
101  do i = 1, num_lagrange
102  if( symmetricmatrixstruc )then
103  indices(countnon0) = np*ndof + i
104  countnon0 = countnon0 + 1
105  else
106  do l = fstrmat%indexL_lagrange(i-1)+1, fstrmat%indexL_lagrange(i)
107  do k = 1, ndof
108  indices(countnon0) = (fstrmat%itemL_lagrange(l)-1)*ndof + k
109  countnon0 = countnon0 + 1
110  enddo
111  enddo
112  endif
113  pointers(np*ndof+i+1) = countnon0
114  enddo
115  endif
116 
118 
119 
123  subroutine set_values_directsolver(hecMAT,fstrMAT)
124 
125  type(hecmwst_matrix) :: hecMAT
126  type (fstrST_matrix_contact_lagrange) :: fstrMAT
127 
128  integer (kind=kint) :: np
129  integer (kind=kint) :: ndof
130  integer (kind=kint) :: num_lagrange
131  integer (kind=kint) :: ierr
132  integer (kind=kint) :: i, j, k, l
133  integer (kind=kint) :: countNon0, locINal, locINd, locINau, locINal_lag, locINau_lag
134 
135  np = hecmat%NP ; ndof = hecmat%NDOF ; num_lagrange = fstrmat%num_lagrange
136 
137  if(allocated(values))deallocate(values)
138  allocate(values(numnon0), stat=ierr)
139  if( ierr /= 0 ) stop " Allocation error, mkl%values "
140  values = 0.0d0
141 
142  countnon0 = 1
143  do i = 1, np
144  do j = 1, ndof
145  if( .not. symmetricmatrixstruc )then
146  do l = hecmat%indexL(i-1)+1, hecmat%indexL(i)
147  do k = 1, ndof
148  locinal = ((l-1)*ndof+j-1)*ndof + k
149  values(countnon0) = hecmat%AL(locinal)
150  countnon0 = countnon0 + 1
151  enddo
152  enddo
153  do k = 1, j-1
154  locind = ((i-1)*ndof+j-1)*ndof + k
155  values(countnon0) = hecmat%D(locind)
156  countnon0 = countnon0 + 1
157  enddo
158  endif
159  do k = j, ndof
160  locind = ((i-1)*ndof+j-1)*ndof + k
161  values(countnon0) = hecmat%D(locind)
162  countnon0 = countnon0 + 1
163  enddo
164  do l = hecmat%indexU(i-1)+1, hecmat%indexU(i)
165  do k = 1, ndof
166  locinau = ((l-1)*ndof+j-1)*ndof + k
167  values(countnon0) = hecmat%AU(locinau)
168  countnon0 = countnon0 + 1
169  enddo
170  enddo
171  if( num_lagrange > 0 )then
172  do l = fstrmat%indexU_lagrange(i-1)+1, fstrmat%indexU_lagrange(i)
173  locinau_lag = (l-1)*ndof + j
174  values(countnon0) = fstrmat%AU_lagrange(locinau_lag)
175  countnon0 = countnon0 + 1
176  enddo
177  endif
178  enddo
179  enddo
180 
181  if( .not.symmetricmatrixstruc .and. num_lagrange > 0 )then
182  do i = 1, num_lagrange
183  do l = fstrmat%indexL_lagrange(i-1)+1, fstrmat%indexL_lagrange(i)
184  do k = 1, ndof
185  locinal_lag = (l-1)*ndof + k
186  values(countnon0) = fstrmat%AL_lagrange(locinal_lag)
187  countnon0 = countnon0 + 1
188  enddo
189  enddo
190  enddo
191  endif
192 
193  end subroutine set_values_directsolver
194 
196  subroutine getapproximateb(ntdf,x,y)
197 
198  integer(kind=kint) :: ntdf
199  integer(kind=kint) :: i, j, k
200  real(kind=kreal) :: x(ntdf)
201  real(kind=kreal) :: y(ntdf)
202 
203  y = 0.0d0
204  do i = 1, ntdf
205  do j = pointers(i), pointers(i+1)-1
206  k = indices(j)
207  y(i) = y(i) + values(j)*x(k)
208  if( symmetricmatrixstruc .and. k/=i )&
209  y(k)=y(k)+values(j)*x(i)
210  enddo
211  enddo
212 
213  end subroutine getapproximateb
214 
215 
216  subroutine checkresidual(hecMAT,fstrMAT)
217  type(hecmwst_matrix) :: hecmat
218  type (fstrst_matrix_contact_lagrange) :: fstrMAT
219  integer(kind=kint) :: ntdf
220  real(kind=kreal), allocatable :: y(:)
221  real(kind=kreal) :: residual_max
222 
223  ntdf = hecmat%NP*hecmat%NDOF + fstrmat%num_lagrange
224 
225  allocate(y(size(hecmat%B)))
226  y = 0.0d0
227  call getapproximateb(ntdf,hecmat%X,y)
228  residual_max=maxval(dabs(y-hecmat%B))
229  write(*,*)' maximum residual = ',residual_max
230  if( dabs(residual_max) >= 1.0d-8) then
231  write(*,*) ' ###Maximum residual exceeded 1.0d-8---Direct Solver### '
232  ! stop
233  endif
234  deallocate(y)
235 
236  end subroutine checkresidual
237 
238 
This module provides functions of reconstructing.
This module defined coomon data and basic structures for analysis.
Definition: m_fstr.f90:15
This module provides functions to set arrays for direct sparse solver \in the case of using standard ...
real(kind=kreal), dimension(:), allocatable values
a
integer(kind=kint), dimension(:), allocatable pointers
ia
subroutine set_pointersandindices_directsolver(hecMAT, fstrMAT, is_sym)
This subroutine sets index arrays for direct sparse solver from those stored \in the matrix structure...
subroutine set_values_directsolver(hecMAT, fstrMAT)
This subroutine sets the array for direct sparse solver that contains \the non-zero items(elements)of...
subroutine getapproximateb(ntdf, x, y)
This subroutine gets the residual vector.
integer(kind=kint), dimension(:), allocatable indices
ja
subroutine checkresidual(hecMAT, fstrMAT)
Structure for Lagrange multiplier-related part of stiffness matrix (Lagrange multiplier-related matri...