FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
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
22contains
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
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 ...
subroutine set_values_directsolver(hecmat, fstrmat)
This subroutine sets the array for direct sparse solver that contains \the non-zero items(elements)of...
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 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...