Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
fractography3d.f90
Go to the documentation of this file.
1 !*********************************************************************
2 !* Illinois Open Source License *
3 !* *
4 !* University of Illinois/NCSA *
5 !* Open Source License *
6 !* *
7 !* Copyright@2008, University of Illinois. All rights reserved. *
8 !* *
9 !* Developed by: *
10 !* *
11 !* Center for Simulation of Advanced Rockets *
12 !* *
13 !* University of Illinois *
14 !* *
15 !* www.csar.uiuc.edu *
16 !* *
17 !* Permission is hereby granted, free of charge, to any person *
18 !* obtaining a copy of this software and associated documentation *
19 !* files (the "Software"), to deal with the Software without *
20 !* restriction, including without limitation the rights to use, *
21 !* copy, modify, merge, publish, distribute, sublicense, and/or *
22 !* sell copies of the Software, and to permit persons to whom the *
23 !* Software is furnished to do so, subject to the following *
24 !* conditions: *
25 !* *
26 !* *
27 !* @ Redistributions of source code must retain the above copyright *
28 !* notice, this list of conditions and the following disclaimers. *
29 !* *
30 !* @ Redistributions in binary form must reproduce the above *
31 !* copyright notice, this list of conditions and the following *
32 !* disclaimers in the documentation and/or other materials *
33 !* provided with the distribution. *
34 !* *
35 !* @ Neither the names of the Center for Simulation of Advanced *
36 !* Rockets, the University of Illinois, nor the names of its *
37 !* contributors may be used to endorse or promote products derived *
38 !* from this Software without specific prior written permission. *
39 !* *
40 !* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, *
41 !* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES *
42 !* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND *
43 !* NONINFRINGEMENT. IN NO EVENT SHALL THE CONTRIBUTORS OR *
44 !* COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
45 !* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, *
46 !* ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
47 !* USE OR OTHER DEALINGS WITH THE SOFTWARE. *
48 !*********************************************************************
49 !* Please acknowledge The University of Illinois Center for *
50 !* Simulation of Advanced Rockets in works and publications *
51 !* resulting from this software or its derivatives. *
52 !*********************************************************************
53 
54 !!****
55 !!
56 !! NAME
57 !! fractography3d
58 !!
59 !! FUNCTION
60 !! This program inputs a model of a solids structure and
61 !! performs an implicit time-dependant analysis of its
62 !! motion based on the Newmark-Wilson scheme.
63 !!
64 !! INPUTS
65 !! File: fractography3d.inp
66 !! File: timestep.in
67 !! Files: <casename>/<casename>.<proc#>.inp
68 !!
69 !! OUTPUTS
70 !! Files: output_<proc#>.dat
71 !!
72 !! USES
73 !! MPI
74 !! readpat
75 !! feminp
76 !! ReadMeshtran
77 !! InitComm1
78 !! InitComm2
79 !! get_mat_stiffness
80 !! implicit_v3d8_mass_consistent
81 !! implicit_v3d8_mass
82 !! implicit_v3d8_me_K
83 !! comp_row_getval
84 !! comp_row_addval
85 !! comp_row_resize
86 !! comp_row_add
87 !! implicit_bc_enforce
88 !! vol_elem_mat
89 !! get_fint
90 !! BS95SETUP
91 !! BS95SOLVE
92 !! BS95FREE
93 !! BS95FINALIZE
94 !! removebcs_meff
95 !! removebcs_pbar
96 !! removebcs_newa
97 !!
98 !!
99 !!****
100 
101 PROGRAM fractography3d
102 
103  USE precision
104 
105  USE comp_row_global
106  USE implicit_global
107 
108  IMPLICIT NONE
109  include 'mpif.h'
110 
111  ! Declare types that you'll need
112  TYPE int_buf
113  INTEGER, DIMENSION(:), POINTER :: rcvbuf
114  END TYPE int_buf
115 
116  ! Declare the variables that you'll need
117  INTEGER :: ios ! error control in opening file
118  INTEGER :: n, nstep ! current timestep & total timesteps
119  INTEGER :: BS95debug = 0 ! debug statements in BlockSolve95
120  REAL(kind=wp) :: delt ! time step size
121  REAL(kind=wp) :: thetaimp, alphaimp, deltaimp ! implicit direct integration constants
122  REAL(kind=wp), ALLOCATABLE, DIMENSION(:) :: fext ! vector containing external forces
123  CHARACTER*180 :: word180 ! dummy string for reading in input files
124  INTEGER :: ifreq ! frequency of output
125  INTEGER :: nnz_k ! number of nonzeros in the K matrix
126  INTEGER :: nnz_m ! number of nonzeros in the M matrix
127  INTEGER :: i, j, m, p, counter ! counters
128  INTEGER :: ngpts = 8 ! number of gauss points per element
129  REAL(kind=wp),ALLOCATABLE,DIMENSION(:) :: tempmg, mg ! Global mass matrix (lumped)
130  REAL(kind=wp),ALLOCATABLE,DIMENSION(:) :: aval_m, aval_k ! M and K matrices in compressed row storage
131  INTEGER,ALLOCATABLE,DIMENSION(:) :: rp_m, cval_m, rp_k, cval_k ! M and K matrices in compressed row storage
132  INTEGER :: nnz_m ! Number of nonzeros in the mass matrix
133  INTEGER :: nstart_km, nrows_km ! Dimensions of parts of M and K matrices
134  REAL(kind=wp),ALLOCATABLE,DIMENSION(:) :: aval_meff ! Effective M matrix in compressed row storage
135  INTEGER,ALLOCATABLE,DIMENSION(:) :: cval_meff, rp_meff ! Effective M matrix in compressed row storage
136  INTEGER :: nnz_meff, nstart_meff, nrows_meff ! Dimensions of parts of effective M matrix
137  REAL(kind=wp),ALLOCATABLE,DIMENSION(:) :: pbar ! Effective load vector
138  REAL(kind=wp),ALLOCATABLE,DIMENSION(:) :: newpbar ! Pbar after removing DOF related to displacement BCs
139  REAL(kind=wp),ALLOCATABLE,DIMENSION(:) :: fint ! Internal force vector
140  REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:) :: S ! Stress used in calculation of internal force vector
141  REAL(kind=wp), DIMENSION(1:3,1:8) :: ri = RESHAPE( &
142  (/-0.577350269189626,-0.577350269189626,-0.577350269189626, &
143  0.577350269189626,-0.577350269189626,-0.577350269189626, &
144  0.577350269189626, 0.577350269189626,-0.577350269189626, &
145  -0.577350269189626, 0.577350269189626,-0.577350269189626, &
146  -0.577350269189626,-0.577350269189626, 0.577350269189626, &
147  0.577350269189626,-0.577350269189626, 0.577350269189626, &
148  0.577350269189626, 0.577350269189626, 0.577350269189626, &
149  -0.577350269189626, 0.577350269189626, 0.577350269189626/),(/3,8/) )
150  REAL(KIND=wp), ALLOCATABLE, DIMENSION(:) :: v ! nodal velocity
151  REAL(KIND=wp), ALLOCATABLE, DIMENSION(:) :: a ! nodal acceleration
152  REAL(kind=wp) :: t ! current time within iteration
153  REAL(kind=wp) :: maxdisp, gmaxdisp ! maximum nodal displacement
154  REAL(kind=wp), ALLOCATABLE, DIMENSION(:) :: ftemp ! temporary vector used to calculate initial acceleration
155  REAL(kind=wp), ALLOCATABLE, DIMENSION(:) :: bufsnd ! MPI send buffer
156  INTEGER, ALLOCATABLE, DIMENSION(:) :: AmountToReceive, AmountToSend ! Amount of stuff that will be received via MPI
157  REAL(kind=wp), ALLOCATABLE, DIMENSION(:) :: finttemp, disptemp ! Temporary arrays used to store info about all NumNp nodes
158  INTEGER :: maxdispnode ! Node that has the maximum displacement
159  REAL(kind=wp) :: tempKval ! temporary value from the K matrix
160  REAL(kind=wp),ALLOCATABLE,DIMENSION(:) :: aval_ktemp ! Temporaray K matrix in compressed row storage
161  INTEGER,ALLOCATABLE,DIMENSION(:) :: cval_ktemp, rp_ktemp ! Temporaray K matrix in compressed row storage
162  INTEGER :: nnz_ktemp ! Number of nonzeros in the temporary K matrix
163  REAL(kind=wp) :: contol ! Convergence tolerance for BS95
164  REAL(kind=wp), DIMENSION(5,9,9) :: dmat ! Material stiffness matrix
165  REAL(kind=wp), DIMENSION(2) :: props ! Material properties needed to find material stiffness matrix
166  REAL(kind=wp), DIMENSION(8) :: xi, eta, zeta ! Gauss point info
167  REAL(kind=wp), DIMENSION(8) :: xiE, etaE, zetaE ! More Gauss point info
168  REAL(kind=wp) :: one = 1.0, three = 3.0
169  REAL(KIND=wp), DIMENSION(1:8,1:9,1:12) :: mixed_map ! Mixed map tensor
170  REAL(KIND=wp), DIMENSION(1:8,1:9,1:9) :: enhanced_map ! Enhanced map tensor
171  INTEGER :: igpt ! Gauss point counter
172  REAL(kind=wp) :: tempval ! value from the K matrix to be sent via MPI
173  INTEGER :: idof, jdof, inode, jnode ! Temporary numbering of nodes during K matrix communication
174  REAL(kind=wp) :: per1 ! Period of the first mode
175  INTEGER :: newnrows_meff ! Number of rows in the Meff matrix after removing BC's
176  INTEGER :: newndim ! Number of DOF after removing DOF associated with imposed displacements
177  INTEGER :: newnstart_meff ! Starting row of the Meff matrix after removing BC's
178  REAL(kind=wp), ALLOCATABLE, DIMENSION(:) :: newa ! Acceleration vector after removing rows pertaining to displacement BCs
179  TYPE(int_buf), ALLOCATABLE, DIMENSION(:) :: Ki, Kj ! i and j indices of nonzeros in K matrix
180  REAL(kind=wp),ALLOCATABLE,DIMENSION(:) :: aval_m_temp ! temporary portion of the M matrix using during construction
181  INTEGER,ALLOCATABLE,DIMENSION(:) :: rp_m_temp, cval_m_temp ! temporary portion of the M matrix using during construction
182  REAL(kind=8) :: elapsedtime
183 
184 
185  ! Initialize MPI
186  CALL mpi_init(ierr)
187  CALL mpi_comm_rank(mpi_comm_world, myid, ierr)
188  CALL mpi_comm_size(mpi_comm_world, nprocs, ierr)
189  WRITE(myid_chr,'(i4.4)') myid
190 
191  ! Start timing of program
192  elapsedtime = -1.0*mpi_wtime()
193 
194 !
195 ! Read geometry, material properties, and FE information
196 !
197 
198  ! Read in timestep information file
199  OPEN(20,file='timestep.in' , form='formatted') !'
200  READ(20,'(a180)') word180
201  READ(20,*) nstep,delt,ifreq,contol ! read in timestep info
202  DO i=1,4
203  READ(20,'(a180)') word180
204  ENDDO
205  READ(20,*) thetaimp,alphaimp,deltaimp ! read in implicit direct integration constants
206  CLOSE(20)
207 
208  !
209  ! -- Open Analysis Deck File
210  scaleload = 1.0
211  OPEN(io_input,file='fractography3d.inp',status='old',iostat=ios)
212  CLOSE(io_input)
213  IF(ios.NE.0)THEN
214  IF(myid.EQ.0) print*, '**WARNING: Unable to find Deck File -- fractography3d.inp'
215  IF(myid.EQ.0) print*, ' - Assuming Patran Neutral File Input Format'
216  CALL readpat()
217 
218  lnumnp = numnp
219  gnumnp = numnp
220  ALLOCATE(local2global(1:numnp))
221  ALLOCATE(global2local(1:numnp))
222  ALLOCATE(nodeproc(1:numnp))
223  DO i = 1, numnp
224  local2global(i) = i
225  global2local(i) = i
226  nodeproc(i) = myid
227  ENDDO
228 
229  ELSE
230  !
231  ! -- Read the control deck file
232  IF(myid.EQ.0) print*,'READING INPUT DECK'
233  CALL feminp()
234  CALL readmeshtran()
235 
236  ENDIF
237 
238 
239 !
240 ! Initialize MPI communications
241 !
242  IF(myid==0) print*,'INITIALIZING COMMUNICATIONS'
243  IF (nprocs > 1) THEN
244  CALL initcomm1()
245  CALL initcomm2()
246  ENDIF
247 
248 !
249 ! Construct the material stiffness matrix for each material
250 !
251  DO i = 1, nummat
252  props(1) = e(i)
253  props(2) = xnu(i)
254  CALL get_mat_stiffness(props,dmat(i,:,:),i)
255  ENDDO
256 
257 
258 !
259 ! Get Gauss point info
260 !
261  xi = (/ one/sqrt(three), -one/sqrt(three), -one/sqrt(three), &
262  one/sqrt(three), one/sqrt(three), -one/sqrt(three), &
263  -one/sqrt(three), one/sqrt(three) /)
264 
265  eta = (/ one/sqrt(three), one/sqrt(three), -one/sqrt(three), &
266  -one/sqrt(three), one/sqrt(three), one/sqrt(three), &
267  -one/sqrt(three), -one/sqrt(three) /)
268 
269  zeta = (/ one/sqrt(three), one/sqrt(three), one/sqrt(three), &
270  one/sqrt(three), -one/sqrt(three), -one/sqrt(three), &
271  -one/sqrt(three), -one/sqrt(three) /)
272 
273  xie = (/ sqrt(three), -sqrt(three), -sqrt(three), &
274  sqrt(three), sqrt(three), -sqrt(three), &
275  -sqrt(three), sqrt(three) /)
276 
277  etae = (/ sqrt(three), sqrt(three), -sqrt(three), &
278  -sqrt(three), sqrt(three), sqrt(three), &
279  -sqrt(three), -sqrt(three) /)
280 
281  zetae = (/ sqrt(three), sqrt(three), sqrt(three), &
282  sqrt(three), -sqrt(three), -sqrt(three), &
283  -sqrt(three), -sqrt(three) /)
284  mixed_map(:,:,:) = 0.0
285  enhanced_map(:,:,:) = 0.0
286  DO igpt = 1, 8
287  mixed_map(igpt,1,1) = eta(igpt)
288  mixed_map(igpt,1,2) = zeta(igpt)
289  mixed_map(igpt,1,3) = eta(igpt) * zeta(igpt)
290  mixed_map(igpt,2,10) = zeta(igpt)
291  mixed_map(igpt,3,12) = eta(igpt)
292  mixed_map(igpt,4,10) = zeta(igpt)
293  mixed_map(igpt,5,4) = xi(igpt)
294  mixed_map(igpt,5,5) = zeta(igpt)
295  mixed_map(igpt,5,6) = xi(igpt) * zeta(igpt)
296  mixed_map(igpt,6,11) = xi(igpt)
297  mixed_map(igpt,7,12) = eta(igpt)
298  mixed_map(igpt,8,11) = xi(igpt)
299  mixed_map(igpt,9,7) = xi(igpt)
300  mixed_map(igpt,9,8) = eta(igpt)
301  mixed_map(igpt,9,9) = xi(igpt) * eta(igpt)
302  enhanced_map(igpt,1,1) = xi(igpt)
303  enhanced_map(igpt,1,2) = xi(igpt) * eta(igpt)
304  enhanced_map(igpt,1,3) = xi(igpt) * zeta(igpt)
305  enhanced_map(igpt,5,4) = eta(igpt)
306  enhanced_map(igpt,5,5) = eta(igpt) * zeta(igpt)
307  enhanced_map(igpt,5,6) = eta(igpt) * xi(igpt)
308  enhanced_map(igpt,9,7) = zeta(igpt)
309  enhanced_map(igpt,9,8) = zeta(igpt) * eta(igpt)
310  enhanced_map(igpt,9,9) = zeta(igpt) * xi(igpt)
311  ENDDO
312  ALLOCATE(aenh(1:9,1:numelv))
313  ALLOCATE(stress(1:ngpts,1:9,1:numelv))
314  aenh(:,:) = 0.0
315  stress(:,:,:) = 0.0
316 
317  ! Set up some constants
318  nstart_km = gnumnp
319  DO m = 1, gnumnp
320  DO i = 1, numnp
321  IF (local2global(i) == m) THEN
322  IF (nodeproc(i) == myid) THEN
323  nstart_km = min(nstart_km,local2global(i))
324  ENDIF
325  ENDIF
326  ENDDO
327  ENDDO
328  nstart_km = 3 * (nstart_km - 1) + 1
329  nrows_km = 3*lnumnp
330 
331 !
332 ! 1. Form the stiffness, and mass matrices
333 !
334 
335 
336  ! Construct the global mass matrix
337  IF(myid==0) print*,'CONSTRUCTING THE MASS MATRIX'
338 
339  ! Construct the mass matrix as a consistent matrix
340  CALL implicit_v3d8_mass_consistent(numelv,numnp,nummat,coor,elconnvol,mattype,ri,rho,elconnvol)
341  nnz_m = nnz_temp
342  ALLOCATE(cval_m(1:nnz_m))
343  ALLOCATE(aval_m(1:nnz_m))
344  ALLOCATE(rp_m(1:3*gnumnp+1))
345  rp_m = rp_temp
346  cval_m = cval_temp
347  aval_m = aval_temp
348  DEALLOCATE(rp_temp)
349  DEALLOCATE(cval_temp)
350  DEALLOCATE(aval_temp)
351 
352  ! Communicate mass to the other procs
353  IF (nprocs > 1) THEN
354  CALL mpi_barrier(mpi_comm_world,ierr)
355  ALLOCATE(req_rcv(1:nprocs))
356  ALLOCATE(req_snd(1:nprocs))
357  ALLOCATE(stat_snd(1:mpi_status_size,1:nprocs))
358  ALLOCATE(stat_rcv(1:mpi_status_size,1:nprocs))
359  DO i = 1, nprocs
360  req_rcv(i) = 0
361  req_snd(i) = 0
362  DO j = 1, mpi_status_size
363  stat_rcv(1,i) = 0
364  stat_snd(j,i) = 0
365  ENDDO
366  ENDDO
367  ALLOCATE(amounttoreceive(1:numcommprocsfrom1))
368  DO i = 1, numcommprocsfrom1
369  CALL mpi_irecv(amounttoreceive(i),1, &
370  mpi_integer,commprocsfrom1(i),10,mpi_comm_world, &
371  req_rcv(i),ierr)
372  ENDDO
373  ALLOCATE(amounttosend(1:numcommprocs1))
374  DO i = 1, numcommprocs1
375  counter = 0
376  DO j = 1, numcommnodes1(i)
377  DO p = 1, 3
378  DO m = 1, 3*gnumnp
379  tempkval = 0.0
380  inode = local2global(commnodes1(i,j))
381  jnode = int((m-0.5)/3)+1
382  idof = p
383  jdof = m - 3*jnode + 3
384  CALL comp_row_getval(3*gnumnp,3*gnumnp,nnz_m,1,rp_m,cval_m,aval_m, &
385  inode*3-3+idof,jnode*3-3+jdof,tempkval)
386  IF (tempkval /= 0.0) THEN
387  counter = counter + 1
388  !print*,myid,Global2Local(inode)*3-3+idof,Global2Local(jnode)*3-3+jdof,tempKval
389  ENDIF
390  ENDDO
391  ENDDO
392  !print*,myid,' sending stuff about node ',Local2Global(CommNodes1(i,j)),' to ',CommProcs1(i),' : ',counter
393  ENDDO
394  amounttosend(i) = counter
395  !CALL MPI_ISEND(AmountToSend(i),1,MPI_INTEGER, &
396  ! CommProcs1(i),10,MPI_COMM_WORLD,req_snd(i),ierr)
397  CALL mpi_send(amounttosend(i),1,mpi_integer, &
398  commprocs1(i),10,mpi_comm_world,ierr)
399  ENDDO
400  !CALL MPI_WAITALL(NumCommProcsFrom1,req_rcv,stat_rcv,ierr)
401  !CALL MPI_WAITALL(NumCommProcs1,req_snd,stat_snd,ierr)
402  CALL mpi_waitall(nprocs,req_rcv,stat_rcv,ierr)
403  !CALL MPI_WAITALL(nprocs,req_snd,stat_snd,ierr)
404 
405  ! Communicate parts of the M matrix to other procs
406  CALL mpi_barrier(mpi_comm_world,ierr)
407  ALLOCATE(frmproc(1:numcommprocsfrom1))
408  DO i = 1, numcommprocsfrom1
409  ALLOCATE(frmproc(i)%rcvbuf(1:3*amounttoreceive(i)))
410  CALL mpi_irecv(frmproc(i)%rcvbuf(1),3*amounttoreceive(i), &
411  mpi_double_precision,commprocsfrom1(i),10,mpi_comm_world, &
412  req_rcv(i),ierr)
413  ENDDO
414  DO i = 1, numcommprocs1
415  ALLOCATE(bufsnd(1:3*amounttosend(i)))
416  bufsnd(:) = 0.0
417  counter = 0
418  DO j = 1, numcommnodes1(i)
419  DO p = 1, 3
420  DO m = 1, 3*gnumnp
421  tempkval = 0.0
422  inode = local2global(commnodes1(i,j))
423  jnode = int((m-0.5)/3)+1
424  idof = p
425  jdof = m - 3*jnode + 3
426  CALL comp_row_getval(3*gnumnp,3*gnumnp,nnz_m,1,rp_m,cval_m,aval_m, &
427  inode*3-3+idof,jnode*3-3+jdof,tempkval)
428  IF (tempkval /= 0.0) THEN
429  counter = counter + 1
430  bufsnd(3*counter-2) = local2global(commnodes1(i,j))*3-3+p
431  bufsnd(3*counter-1) = m
432  bufsnd(3*counter) = tempkval
433  ENDIF
434  ENDDO
435  ENDDO
436  ENDDO
437  !print*,myid,' sending to ',CommProcs1(i),' : ',bufsnd(:)
438  !CALL MPI_ISEND(bufsnd,3*AmountToSend(i),MPI_DOUBLE_PRECISION, &
439  ! CommProcs1(i),10,MPI_COMM_WORLD,req_snd(i),ierr)
440  CALL mpi_send(bufsnd,3*amounttosend(i),mpi_double_precision, &
441  commprocs1(i),10,mpi_comm_world,ierr)
442  DEALLOCATE(bufsnd)
443  ENDDO
444  !CALL MPI_WAITALL(NumCommProcsFrom1,req_rcv,stat_rcv,ierr)
445  !CALL MPI_WAITALL(NumCommProcs1,req_snd,stat_snd,ierr)
446  CALL mpi_waitall(nprocs,req_rcv,stat_rcv,ierr)
447  !CALL MPI_WAITALL(nprocs,req_snd,stat_snd,ierr)
448  DEALLOCATE(req_rcv)
449  DEALLOCATE(req_snd)
450  DEALLOCATE(stat_snd)
451  DEALLOCATE(stat_rcv)
452  ENDIF
453 
454  IF (nprocs > 1) THEN
455  DO i = 1, numcommprocsfrom1
456  !print*,myid,' received from ',CommProcsFrom1(i), ' : ',frmproc(i)%rcvbuf(:)
457  DO j = 1, amounttoreceive(i)
458  CALL comp_row_addval(3*gnumnp,3*gnumnp,nnz_m,1,rp_m,cval_m,aval_m, &
459  int(frmproc(i)%rcvbuf(3*j-2)),int(frmproc(i)%rcvbuf(3*j-1)),frmproc(i)%rcvbuf(3*j))
460  IF (nnz_temp /= nnz_m) THEN
461  DEALLOCATE(cval_m,aval_m)
462  ALLOCATE(cval_m(nnz_temp),aval_m(nnz_temp))
463  nnz_m = nnz_temp
464  rp_m = rp_temp
465  cval_m = cval_temp
466  ENDIF
467  aval_m = aval_temp
468  DEALLOCATE(rp_temp,cval_temp,aval_temp)
469  ENDDO
470  ENDDO
471  DEALLOCATE(frmproc)
472  DEALLOCATE(amounttoreceive)
473  DEALLOCATE(amounttosend)
474  ENDIF
475 
476  ! Resize the matrix to the size needed on this proc
477  CALL comp_row_resize(3*gnumnp,3*gnumnp,nnz_m,1,rp_m,cval_m,aval_m,nstart_km,nrows_km)
478  DEALLOCATE(rp_m,cval_m,aval_m)
479  ALLOCATE(rp_m(1:nrows_km+1))
480  ALLOCATE(cval_m(1:nnz_temp))
481  ALLOCATE(aval_m(1:nnz_temp))
482  nnz_m = nnz_temp
483  rp_m = rp_temp
484  cval_m = cval_temp
485  aval_m = aval_temp
486  DEALLOCATE(rp_temp,cval_temp,aval_temp)
487 
488  ! Construct the mass matrix as lumped
489  ALLOCATE(tempmg(1:numnp))
490  ALLOCATE(mg(1:3*lnumnp))
491  ALLOCATE(cval_m_temp(1:3*lnumnp))
492  ALLOCATE(aval_m_temp(1:3*lnumnp))
493  ALLOCATE(rp_m_temp(1:3*lnumnp+1))
494  tempmg(1:numnp) = 0.0
495  CALL implicit_v3d8_mass(numelv,numnp,nummat,coor,elconnvol,mattype,ri,rho,tempmg,1,numelv)
496  DO i = 1, numnp
497  tempmg(i) = 1.0 / tempmg(i) ! Get the mass matrix, not its inverse
498  ENDDO
499 
500  IF (nprocs > 1) THEN
501  CALL mpi_barrier(mpi_comm_world,ierr)
502  !ALLOCATE(req_rcv(1:NumCommProcsFrom1))
503  !ALLOCATE(req_snd(1:NumCommProcs1))
504  !ALLOCATE(stat_snd(1:MPI_STATUS_SIZE,1:NumCommProcs1))
505  !ALLOCATE(stat_rcv(1:MPI_STATUS_SIZE,1:NumCommProcsFrom1))
506  ALLOCATE(req_rcv(1:nprocs))
507  ALLOCATE(req_snd(1:nprocs))
508  ALLOCATE(stat_snd(1:mpi_status_size,1:nprocs))
509  ALLOCATE(stat_rcv(1:mpi_status_size,1:nprocs))
510  DO i = 1, nprocs
511  req_rcv(i) = 0
512  req_snd(i) = 0
513  DO j = 1, mpi_status_size
514  stat_rcv(1,i) = 0
515  stat_snd(j,i) = 0
516  ENDDO
517  ENDDO
518  ALLOCATE(frmproc(1:numcommprocsfrom1))
519  DO i = 1, numcommprocsfrom1
520  ALLOCATE(frmproc(i)%rcvbuf(1:2*numcommnodesfrom1(i)))
521  CALL mpi_irecv(frmproc(i)%rcvbuf(1),2*numcommnodesfrom1(i), &
522  mpi_double_precision,commprocsfrom1(i),10,mpi_comm_world, &
523  req_rcv(i),ierr)
524  ENDDO
525  DO i = 1, numcommprocs1
526  ALLOCATE(bufsnd(2*numcommnodes1(i)))
527  DO j = 1, numcommnodes1(i)
528  bufsnd(2*j-1) = local2global(commnodes1(i,j))
529  bufsnd(2*j) = tempmg(commnodes1(i,j))
530  ENDDO
531  !CALL MPI_ISEND(bufsnd,2*NumCommNodes1(i),MPI_DOUBLE_PRECISION, &
532  ! CommProcs1(i),10,MPI_COMM_WORLD,req_snd(i),ierr)
533  CALL mpi_send(bufsnd,2*numcommnodes1(i),mpi_double_precision, &
534  commprocs1(i),10,mpi_comm_world,ierr)
535  DEALLOCATE(bufsnd)
536  ENDDO
537  !CALL MPI_WAITALL(NumCommProcsFrom1,req_rcv,stat_rcv,ierr)
538  !CALL MPI_WAITALL(NumCommProcs1,req_snd,stat_snd,ierr)
539  CALL mpi_waitall(nprocs,req_rcv,stat_rcv,ierr)
540  !CALL MPI_WAITALL(nprocs,req_snd,stat_snd,ierr)
541  DEALLOCATE(req_rcv)
542  DEALLOCATE(req_snd)
543  DEALLOCATE(stat_snd)
544  DEALLOCATE(stat_rcv)
545 
546  ! Add the communicated mass to the current mass
547  DO i = 1, numcommprocsfrom1
548  DO j = 1, numcommnodesfrom1(i)
549  tempmg(global2local(int(frmproc(i)%rcvbuf(2*j-1)))) = &
550  tempmg(global2local(int(frmproc(i)%rcvbuf(2*j-1)))) + frmproc(i)%rcvbuf(2*j)
551  ENDDO
552  ENDDO
553 
554  DEALLOCATE(frmproc)
555  ENDIF
556 
557  ! Put the mass into compressed row format
558  counter = 0
559  DO m = 1, gnumnp
560  DO i = 1, numnp
561  IF (local2global(i) == m) THEN
562  IF (nodeproc(i) == myid) THEN
563  DO j = 1, 3
564  counter = counter + 1
565  mg(counter) = tempmg(i)
566  ENDDO
567  ENDIF
568  ENDIF
569  ENDDO
570  ENDDO
571  DO i = 1, 3*lnumnp
572  rp_m_temp(i) = i-1
573  cval_m_temp(i) = i-1 + (nstart_km - 1)
574  aval_m_temp(i) = mg(i)
575  ENDDO
576  rp_m_temp(3*lnumnp+1) = 3*lnumnp
577  DEALLOCATE(tempmg)
578  DEALLOCATE(mg)
579 
580  ! Average the consistent and lumped mass matrices to get a higher-order mass matrix
581  CALL comp_row_add(3*lnumnp,3*gnumnp,nrows_km,nrows_km,nnz_m,3*lnumnp,1, &
582  rp_m,cval_m,0.5*aval_m,1,rp_m_temp,cval_m_temp,0.5*aval_m_temp)
583  DEALLOCATE(rp_m,cval_m,aval_m)
584  DEALLOCATE(rp_m_temp,cval_m_temp,aval_m_temp)
585  ALLOCATE(rp_m(nrows_km+1),cval_m(nnz_temp),aval_m(nnz_temp))
586  nnz_m = nnz_temp
587  rp_m = rp_temp
588  cval_m = cval_temp
589  aval_m = aval_temp
590  DEALLOCATE(rp_temp,cval_temp,aval_temp)
591 
592  !! Output the matrix to matlab files for analysis
593  !OPEN(30,FILE='m_'//myid_chr//'.m',FORM='formatted')
594  !counter = 0
595  !DO i = 1, nrows_km
596  ! DO j = rp_m(i)+1, rp_m(i+1)
597  ! counter = counter + 1
598  ! m = cval_m(counter)+1
599  ! WRITE(30,*) 'A(', i+nstart_km-1, ',', m, ') = ',aval_m(counter), ';'
600  ! ENDDO
601  !ENDDO
602  !CLOSE(30)
603 
604 
605  ! Construct the global stiffness matrix
606  IF(myid==0) print*,'CONSTRUCTING THE STIFFNESS MATRIX'
607  CALL implicit_v3d8_me_k(coor, elconnvol, dmat, numnp, 1, numelv, numelv, &
608  mattype, nummat, enhanced_map, mixed_map)
609  ALLOCATE(rp_k(1:3*gnumnp+1))
610  ALLOCATE(cval_k(1:nnz_temp))
611  ALLOCATE(aval_k(1:nnz_temp))
612  nnz_k = nnz_temp
613  rp_k = rp_temp
614  cval_k = cval_temp
615  aval_k = aval_temp
616  DEALLOCATE(rp_temp,cval_temp,aval_temp)
617 
618  ! Communicate the K matrix to other procs
619  IF (nprocs > 1) THEN
620  CALL mpi_barrier(mpi_comm_world,ierr)
621  ALLOCATE(req_rcv(1:nprocs))
622  ALLOCATE(req_snd(1:nprocs))
623  ALLOCATE(stat_snd(1:mpi_status_size,1:nprocs))
624  ALLOCATE(stat_rcv(1:mpi_status_size,1:nprocs))
625  DO i = 1, nprocs
626  req_rcv(i) = 0
627  req_snd(i) = 0
628  DO j = 1, mpi_status_size
629  stat_rcv(1,i) = 0
630  stat_snd(j,i) = 0
631  ENDDO
632  ENDDO
633  ALLOCATE(amounttoreceive(1:numcommprocsfrom1))
634  DO i = 1, numcommprocsfrom1
635  CALL mpi_irecv(amounttoreceive(i),1, &
636  mpi_integer,commprocsfrom1(i),10,mpi_comm_world, &
637  req_rcv(i),ierr)
638  ENDDO
639  ALLOCATE(amounttosend(1:numcommprocs1))
640  DO i = 1, numcommprocs1
641  counter = 0
642  DO j = 1, numcommnodes1(i)
643  DO p = 1, 3
644  DO m = 1, 3*gnumnp
645  tempkval = 0.0
646  inode = local2global(commnodes1(i,j))
647  jnode = int((m-0.5)/3)+1
648  idof = p
649  jdof = m - 3*jnode + 3
650  CALL comp_row_getval(3*gnumnp,3*gnumnp,nnz_k,1,rp_k,cval_k,aval_k, &
651  inode*3-3+idof,jnode*3-3+jdof,tempkval)
652  IF (tempkval /= 0.0) THEN
653  counter = counter + 1
654  !print*,myid,Global2Local(inode)*3-3+idof,Global2Local(jnode)*3-3+jdof,tempKval
655  ENDIF
656  ENDDO
657  ENDDO
658  !print*,myid,' sending stuff about node ',Local2Global(CommNodes1(i,j)), &
659  !' to ',CommProcs1(i),' : ',counter
660  ENDDO
661  amounttosend(i) = counter
662  !CALL MPI_ISEND(AmountToSend(i),1,MPI_INTEGER, &
663  ! CommProcs1(i),10,MPI_COMM_WORLD,req_snd(i),ierr)
664  CALL mpi_send(amounttosend(i),1,mpi_integer, &
665  commprocs1(i),10,mpi_comm_world,ierr)
666  ENDDO
667  !CALL MPI_WAITALL(NumCommProcsFrom1,req_rcv,stat_rcv,ierr)
668  !CALL MPI_WAITALL(NumCommProcs1,req_snd,stat_snd,ierr)
669  CALL mpi_waitall(nprocs,req_rcv,stat_rcv,ierr)
670  !CALL MPI_WAITALL(nprocs,req_snd,stat_snd,ierr)
671 
672 
673  ! Communicate i locations of nonzeros in K matrix to other procs
674  CALL mpi_barrier(mpi_comm_world,ierr)
675  ALLOCATE(ki(1:numcommprocsfrom1))
676  DO i = 1, numcommprocsfrom1
677  ALLOCATE(ki(i)%rcvbuf(1:amounttoreceive(i)))
678  CALL mpi_irecv(ki(i)%rcvbuf(1),amounttoreceive(i), &
679  mpi_integer,commprocsfrom1(i),10,mpi_comm_world, &
680  req_rcv(i),ierr)
681  ENDDO
682  DO i = 1, numcommprocs1
683  ALLOCATE(bufsnd(1:amounttosend(i)))
684  counter = 0
685  DO j = 1, numcommnodes1(i)
686  DO p = 1, 3
687  DO m = 1, 3*gnumnp
688  tempkval = 0.0
689  inode = local2global(commnodes1(i,j))
690  jnode = int((m-0.5)/3)+1
691  idof = p
692  jdof = m - 3*jnode + 3
693  CALL comp_row_getval(3*gnumnp,3*gnumnp,nnz_k,1,rp_k,cval_k,aval_k, &
694  inode*3-3+idof,jnode*3-3+jdof,tempkval)
695  IF (tempkval /= 0.0) THEN
696  counter = counter + 1
697  bufsnd(counter) = local2global(commnodes1(i,j))*3-3+p
698  ENDIF
699  ENDDO
700  ENDDO
701  ENDDO
702  !CALL MPI_ISEND(INT(bufsnd(:)),AmountToSend(i),MPI_INTEGER, &
703  ! CommProcs1(i),10,MPI_COMM_WORLD,req_snd(i),ierr)
704  CALL mpi_send(int(bufsnd(:)),amounttosend(i),mpi_integer, &
705  commprocs1(i),10,mpi_comm_world,ierr)
706  DEALLOCATE(bufsnd)
707  ENDDO
708  CALL mpi_waitall(nprocs,req_rcv,stat_rcv,ierr)
709  !CALL MPI_WAITALL(nprocs,req_snd,stat_snd,ierr)
710 
711  ! Communicate j locations of nonzeros in K matrix to other procs
712  CALL mpi_barrier(mpi_comm_world,ierr)
713  ALLOCATE(kj(1:numcommprocsfrom1))
714  DO i = 1, numcommprocsfrom1
715  ALLOCATE(kj(i)%rcvbuf(1:amounttoreceive(i)))
716  CALL mpi_irecv(kj(i)%rcvbuf(1),amounttoreceive(i), &
717  mpi_integer,commprocsfrom1(i),10,mpi_comm_world, &
718  req_rcv(i),ierr)
719  ENDDO
720  DO i = 1, numcommprocs1
721  ALLOCATE(bufsnd(1:amounttosend(i)))
722  counter = 0
723  DO j = 1, numcommnodes1(i)
724  DO p = 1, 3
725  DO m = 1, 3*gnumnp
726  tempkval = 0.0
727  inode = local2global(commnodes1(i,j))
728  jnode = int((m-0.5)/3)+1
729  idof = p
730  jdof = m - 3*jnode + 3
731  CALL comp_row_getval(3*gnumnp,3*gnumnp,nnz_k,1,rp_k,cval_k,aval_k, &
732  inode*3-3+idof,jnode*3-3+jdof,tempkval)
733  IF (tempkval /= 0.0) THEN
734  counter = counter + 1
735  bufsnd(counter) = m
736  ENDIF
737  ENDDO
738  ENDDO
739  ENDDO
740  !CALL MPI_ISEND(INT(bufsnd(:)),AmountToSend(i),MPI_INTEGER, &
741  ! CommProcs1(i),10,MPI_COMM_WORLD,req_snd(i),ierr)
742  CALL mpi_send(int(bufsnd(:)),amounttosend(i),mpi_integer, &
743  commprocs1(i),10,mpi_comm_world,ierr)
744  DEALLOCATE(bufsnd)
745  ENDDO
746  CALL mpi_waitall(nprocs,req_rcv,stat_rcv,ierr)
747  !CALL MPI_WAITALL(nprocs,req_snd,stat_snd,ierr)
748 
749  ! Communicate values of nonzeros in K matrix to other procs
750  CALL mpi_barrier(mpi_comm_world,ierr)
751  ALLOCATE(frmproc(1:numcommprocsfrom1))
752  DO i = 1, numcommprocsfrom1
753  ALLOCATE(frmproc(i)%rcvbuf(1:amounttoreceive(i)))
754  CALL mpi_irecv(frmproc(i)%rcvbuf(1),amounttoreceive(i), &
755  mpi_double_precision,commprocsfrom1(i),10,mpi_comm_world, &
756  req_rcv(i),ierr)
757  ENDDO
758  DO i = 1, numcommprocs1
759  ALLOCATE(bufsnd(1:amounttosend(i)))
760  counter = 0
761  DO j = 1, numcommnodes1(i)
762  DO p = 1, 3
763  DO m = 1, 3*gnumnp
764  tempkval = 0.0
765  inode = local2global(commnodes1(i,j))
766  jnode = int((m-0.5)/3)+1
767  idof = p
768  jdof = m - 3*jnode + 3
769  CALL comp_row_getval(3*gnumnp,3*gnumnp,nnz_k,1,rp_k,cval_k,aval_k, &
770  inode*3-3+idof,jnode*3-3+jdof,tempkval)
771  IF (tempkval /= 0.0) THEN
772  counter = counter + 1
773  bufsnd(counter) = tempkval
774  ENDIF
775  ENDDO
776  ENDDO
777  ENDDO
778  !CALL MPI_ISEND(bufsnd,AmountToSend(i),MPI_DOUBLE_PRECISION, &
779  ! CommProcs1(i),10,MPI_COMM_WORLD,req_snd(i),ierr)
780  CALL mpi_send(bufsnd,amounttosend(i),mpi_double_precision, &
781  commprocs1(i),10,mpi_comm_world,ierr)
782  DEALLOCATE(bufsnd)
783  ENDDO
784  CALL mpi_waitall(nprocs,req_rcv,stat_rcv,ierr)
785  !CALL MPI_WAITALL(nprocs,req_snd,stat_snd,ierr)
786 
787  DEALLOCATE(req_rcv)
788  DEALLOCATE(req_snd)
789  DEALLOCATE(stat_snd)
790  DEALLOCATE(stat_rcv)
791 
792  ENDIF
793 
794  CALL comp_row_resize(3*gnumnp,3*gnumnp,nnz_k,1,rp_k,cval_k,aval_k,nstart_km,nrows_km)
795  DEALLOCATE(rp_k,cval_k,aval_k)
796  ALLOCATE(rp_k(1:nrows_km+1))
797  ALLOCATE(cval_k(1:nnz_temp))
798  ALLOCATE(aval_k(1:nnz_temp))
799  nnz_k = nnz_temp
800  rp_k = rp_temp
801  cval_k = cval_temp
802  aval_k = aval_temp
803  DEALLOCATE(rp_temp,cval_temp,aval_temp)
804 
805 
806  IF (nprocs > 1) THEN
807  DO i = 1, numcommprocsfrom1
808  !print*,myid,' received from ',CommProcsFrom1(i), ' : ',frmproc(i)%rcvbuf(:)
809  DO j = 1, amounttoreceive(i)
810  CALL comp_row_addval(3*gnumnp,nrows_km,nnz_k,nstart_km,rp_k,cval_k,aval_k, &
811  ki(i)%rcvbuf(j),kj(i)%rcvbuf(j),frmproc(i)%rcvbuf(j))
812  IF (nnz_temp /= nnz_k) THEN
813  DEALLOCATE(cval_k,aval_k)
814  ALLOCATE(cval_k(nnz_temp),aval_k(nnz_temp))
815  nnz_k = nnz_temp
816  rp_k = rp_temp
817  ENDIF
818  cval_k = cval_temp
819  aval_k = aval_temp
820  DEALLOCATE(rp_temp,cval_temp,aval_temp)
821  ENDDO
822  ENDDO
823  DEALLOCATE(frmproc)
824  DEALLOCATE(ki)
825  DEALLOCATE(kj)
826  ENDIF
827 
828 
829 ! .
830 ! 2. Initialize displacement and velocity array: Uo,Uo
831 !
832  IF(myid==0) print*,'INITIALIZING DISPLACEMENT AND VELOCITY'
833  ALLOCATE(disp(1:3*lnumnp),v(1:3*lnumnp),a(1:3*lnumnp))
834  disp(1:3*lnumnp) = 0.0
835  v(1:3*lnumnp) = 0.0
836  IF (ALLOCATED(node_flag)) THEN
837  ELSE
838  ALLOCATE(node_flag(1:numnp,1:3))
839  ALLOCATE(boundary_value(1:numnp,1:3))
840  node_flag(1:numnp,1:3) = 0
841  DO i = 1, numbound
842  DO j = 1, 3
843 !!$ IF (BC_ID(j+1,i) == 0) THEN ! displacement BC
844 !!$ node_flag(BC_ID(1,i),j) = 8
845 !!$ boundary_value(BC_ID(1,i),j) = BC_VAL(j,i)
846  IF ( global%BCFlag(j+1,i) == 0 ) THEN
847  node_flag(global%BCFlag(1,i),j) = 8
848  boundary_value(global%BCFlag(1,i),j) = global%BCvalue(j,i)
849  ENDIF
850 !!$ IF (BC_ID(j+1,i) == 1) THEN ! force BC
851 !!$ node_flag(BC_ID(1,i),j) = 7
852 !!$ boundary_value(BC_ID(1,i),j) = BC_VAL(j,i)
853  IF ( global%BCFlag(j+1,i) == 1 ) THEN
854  node_flag(global%BCFlag(1,i),j) = 7
855  boundary_value(global%BCFlag(1,i),j) = global%BCvalue(j,i)
856  ENDIF
857  ENDDO
858  ENDDO
859  ENDIF
860  CALL implicit_bc_enforce(numnp,lnumnp,disp,a,v,node_flag,boundary_value,0,myid)
861 
862 
863  !
864  ! Form the material compliance matrix
865  !
866 
867  IF(ALLOCATED(c)) DEALLOCATE(c)
868  IF(ALLOCATED(ci)) DEALLOCATE(ci)
869  ALLOCATE(c(1:9,1:nummat),ci(1:9,1:nummat))
870  c(1:9,1:nummat) = 0.0
871  ci(1:9,1:nummat) = 0.0
872  CALL vol_elem_mat(e,xnu,c,ci,nummat)
873 
874 
875 ! ..
876 ! 3. Calculate the initial acceleration array: Uo
877 !
878 
879 
880  ! Use applied forces and mass matrix
881  ! i.e. solve {F} - [K]*{disp} = [M]*{a} for {a}
882 
883  IF(myid==0) print*,'CALCULATING INITIAL ACCELERATION'
884 
885  ALLOCATE(fint(1:3*lnumnp))
886  CALL get_fint(3*gnumnp,nrows_km,nnz_k,nstart_km,rp_k,cval_k,aval_k,disp,fint)
887 
888  ! Construct the external force vector
889  ALLOCATE(fext(1:3*lnumnp))
890  fext(1:3*lnumnp) = 0.0
891  counter = 0
892  DO m = 1, gnumnp
893  DO i = 1, numnp
894  IF (local2global(i) == m) THEN
895  DO j = 1, 3
896  IF (nodeproc(i) == myid) THEN
897  counter = counter + 1
898  IF(node_flag(i,j) == 7) THEN ! Imposed force
899  fext(counter) = boundary_value(i,j)
900  !print*,myid,' external force on global node ',Local2Global(i),' in direction ',j,' of value ',fext(counter)
901  ENDIF
902  ENDIF
903  ENDDO
904  ENDIF
905  ENDDO
906  ENDDO
907  ALLOCATE(ftemp(1:3*lnumnp))
908  ftemp(1:3*lnumnp) = 0.0
909 
910 
911  ! Construct the effective force vector
912  DO i = 1, 3*lnumnp
913  ftemp(i) = fext(i) - fint(i)
914  !print*,myid,i,ftemp(i)
915  ENDDO
916  DEALLOCATE(fint)
917 
918  ! Remove rows that are constrained by displacement BC's
919  counter = 0
920  DO m = 1, gnumnp
921  DO i = 1, numnp
922  IF (local2global(i) == m) THEN
923  IF (nodeproc(i)==myid) THEN
924  DO j = 1, 3
925  counter = counter + 1
926  IF(node_flag(i,j) == 8) THEN ! Imposed constant nodal displacement
927  DO n = rp_m(counter)+1,rp_m(counter+1)
928  IF (cval_m(n)+1 /= counter+nstart_km-1) THEN
929  aval_m(n) = 0.0
930  ELSE
931  aval_m(n) = 1.0
932  ENDIF
933  ENDDO
934  ENDIF
935  ENDDO
936  ENDIF
937  ENDIF
938  ENDDO
939  ENDDO
940  counter = 0
941  DO m = 1, gnumnp
942  DO i = 1, numnp
943  IF (local2global(i) == m) THEN
944  IF (nodeproc(i)==myid) THEN
945  DO j = 1, 3
946  counter = counter + 1
947  IF(node_flag(i,j) == 8) THEN ! Imposed constant nodal displacement
948  ftemp(counter) = 0.0 ! If constant displacement, acceleration = 0
949  ENDIF
950  ENDDO
951  ENDIF
952  ENDIF
953  ENDDO
954  ENDDO
955 
956 
957 !
958 ! 3.1 Solve for acceleration array
959 !
960  ! Solve using BlockSolve95 functions
961  a(:) = 0.0 ! Initial guess
962  CALL bs95setup(3*gnumnp,nnz_m,nstart_km-1,nrows_km,rp_m,cval_m,aval_m,1,bs95debug)
963  CALL bs95solve(3*lnumnp,ftemp,a,contol,bs95debug)
964  CALL bs95free(bs95debug)
965  DEALLOCATE(ftemp)
966 
967 
968 !
969 ! Make sure initial conditions conform to BC's
970 !
971  CALL implicit_bc_enforce(numnp,lnumnp,disp,a,v,node_flag,boundary_value,0,myid)
972 
973 
974 ! _ _
975 ! Form effective mass matrix M: M = M + gamma*Dt*C + beta*Dt^2*K
976 !
977 
978  IF(myid==0) print*,'FORMING EFFECTIVE MASS MATRIX'
979  CALL comp_row_add(3*lnumnp,3*gnumnp,nrows_km,nrows_km,nnz_k,nnz_m,1, &
980  rp_k,cval_k,alphaimp*delt*delt*aval_k,1,rp_m,cval_m,aval_m)
981  ALLOCATE(rp_meff(nrows_km+1),cval_meff(nnz_temp),aval_meff(nnz_temp))
982  nnz_meff = nnz_temp
983  nstart_meff = nstart_km
984  nrows_meff = nrows_km
985  rp_meff = rp_temp
986  cval_meff = cval_temp
987  aval_meff = aval_temp
988  DEALLOCATE(rp_temp,cval_temp,aval_temp)
989 
990  ! Remove rows that are constrained by displacement BC's
991  CALL removebcs_meff(3*gnumnp,nrows_meff,nnz_meff,nstart_meff,rp_meff,cval_meff,aval_meff, &
992  newnrows_meff,newnstart_meff,newndim)
993  DEALLOCATE(rp_meff,cval_meff,aval_meff)
994  nrows_meff = newnrows_meff
995  nstart_meff = newnstart_meff
996  ALLOCATE(rp_meff(nrows_meff+1),cval_meff(nnz_temp),aval_meff(nnz_temp))
997  nnz_meff = nnz_temp
998  rp_meff = rp_temp
999  cval_meff = cval_temp
1000  aval_meff = aval_temp
1001  DEALLOCATE(rp_temp,cval_temp,aval_temp)
1002 
1003  !! Output matrix to matlab files for analysis
1004  !OPEN(30,FILE='meff_'//myid_chr//'.m',FORM='formatted')
1005  !counter = 0
1006  !DO i = 1, nrows_meff
1007  ! DO j = rp_meff(i)+1, rp_meff(i+1)
1008  ! counter = counter + 1
1009  ! m = cval_meff(counter)+1
1010  ! WRITE(30,*) 'A(', i+nstart_meff-1, ',', m, ') = ',aval_meff(counter), ';'
1011  ! ENDDO
1012  !ENDDO
1013  !CLOSE(30)
1014  !OPEN(30,FILE='k_'//myid_chr//'.m',FORM='formatted')
1015  !counter = 0
1016  !DO i = 1, nrows_km
1017  ! DO j = rp_k(i)+1, rp_k(i+1)
1018  ! counter = counter + 1
1019  ! m = cval_k(counter)+1
1020  ! WRITE(30,*) 'A(', i+nstart_km-1, ',', m, ') = ',aval_k(counter), ';'
1021  ! ENDDO
1022  !ENDDO
1023  !CLOSE(30)
1024 
1025  CALL bs95setup(newndim,nnz_meff,nstart_meff-1,nrows_meff,rp_meff,cval_meff,aval_meff,1,bs95debug)
1026 
1027 
1028 !
1029 ! Allocate variables for use within the time loop
1030 !
1031  ALLOCATE(pbar(1:3*lnumnp))
1032  ALLOCATE(fint(1:3*lnumnp))
1033  ALLOCATE(s(1:6,1:ngpts,1:numelv))
1034  ALLOCATE(newpbar(1:nrows_meff))
1035  ALLOCATE(newa(1:nrows_meff))
1036  newa(:) = 0.0 ! Initial guess
1037 
1038 
1039 !
1040 ! Open output files
1041 !
1042  OPEN(30,file='output_'//myid_chr//'.dat',form='formatted')
1043 
1044  !DO i = 1, NumNp
1045  ! IF(NodeProc(i) == myid) THEN
1046  ! counter = counter + 1
1047  ! !IF(Local2Global(i) == 41)THEN
1048  ! !IF(Local2Global(i) == 88) THEN
1049  ! IF(Local2Global(i) == 718) THEN
1050  ! OPEN(40,FILE='waveprop_'//myid_chr//'.dat',FORM='formatted')
1051  ! EXIT
1052  ! ENDIF
1053  ! ENDIF
1054  !ENDDO
1055 
1056 
1057 
1058 !
1059 ! For each time step
1060 
1061  IF (myid==0) THEN
1062  print*,'BEGINNING TIMESTEPPING'
1063  print*,' '
1064  print*,' step time maxdisp'
1065  print*,'------------------------------------'
1066  ENDIF
1067  t = 0.0
1068  DO n=1,nstep
1069 
1070  t = t + delt
1071  !IF(myid==0) print*,'timestep ',n,' of ',nstep,' delt = ',delt,' t = ',t
1072 
1073 
1074 ! Calculate displacement and velocity predictors:
1075 ! ~ . ..
1076 ! x = x + Dt * x + 1/2 * Dt^2 * (1-2*beta) * x
1077 ! k+1 k k k
1078 !
1079 ! ~. . ..
1080 ! x = x + Dt * (1-gamma) * x
1081 ! k+1 k k
1082 !
1083 
1084  DO i = 1, 3*lnumnp
1085  disp(i) = disp(i) + delt * v(i) + 0.5*delt*delt * (1.0 - 2.0*alphaimp) * a(i)
1086  v(i) = v(i) + delt * (1.0 - deltaimp) * a(i)
1087  ENDDO
1088 
1089 
1090 ! Calculate effective loads at time t + Dt:
1091 !
1092 ! _ .
1093 ! P = P - C x - I
1094 ! k+1 k+1 k
1095 !
1096 ! where I is the internal force vector
1097 !
1098  ! Get the internal force vector
1099  CALL get_fint(3*gnumnp,nrows_km,nnz_k,nstart_km,rp_k,cval_k,aval_k,disp,fint)
1100 
1101  ! Construct the load vector
1102  ! BEGIN TIME-DEPENDENT FORCE CODE MODIFICATIONS
1103  !per1 = 359.159245157598
1104  DO i = 1,3*lnumnp
1105  !if(t < per1) then
1106  ! pbar(i) = fext(i)*SIN(3.14159*t/per1) - fint(i)
1107  !else
1108  ! pbar(i) = -fint(i)
1109  !endif
1110  pbar(i) = fext(i) - fint(i)
1111  ENDDO
1112  ! END TIME-DEPENDENT FORCE CODE MODIFICATIONS
1113 
1114  ! Remove rows that are constrained by displacement BC's
1115  CALL removebcs_pbar(nstart_km,3*lnumnp,pbar,nrows_meff,newpbar)
1116 
1117 
1118 !
1119 ! Solve for acceleration at time t + Dt
1120 ! _ _
1121 ! M a = P
1122 ! k+1 k+1
1123 !
1124  ! Solve using BlockSolve95 functions
1125  !newa(:) = 0.0
1126  CALL bs95solve(nrows_meff,newpbar,newa,contol,bs95debug)
1127 !!$ CALL BS95SOLVE(nrows_meff,pbar,a,contol,BS95debug)
1128 
1129 
1130  ! Put newa back into a
1131  CALL removebcs_newa(nstart_km,3*lnumnp,a,nrows_meff,newa)
1132 
1133 
1134 ! Calculate displacement and velocity correctors:
1135 ! ~ ..
1136 ! x = x + beta * Dt^2 * x
1137 ! k+1 k+1 k+1
1138 !
1139 ! . ~. ..
1140 ! x = x + gamma * Dt * x
1141 ! k+1 k+1 k+1
1142 !
1143  DO i = 1, 3*lnumnp
1144  disp(i) = disp(i) + alphaimp*delt*delt*a(i)
1145  v(i) = v(i) + deltaimp*delt*a(i)
1146  ENDDO
1147 
1148 
1149 !
1150 ! 4. Apply boundary conditions and output results
1151 !
1152  CALL implicit_bc_enforce(numnp,lnumnp,disp,v,a,node_flag,boundary_value,t,myid)
1153 
1154 
1155 !
1156 ! Compute maximum nodal displacement
1157 !
1158  maxdispnode = 0
1159  maxdisp = 0.0
1160  i = 0
1161  DO j = 1, numnp
1162  IF (nodeproc(j) == myid) THEN
1163  i = i + 1
1164  IF ( maxdisp < sqrt(disp(3*(i-1)+1)**2.0 + disp(3*(i-1)+2)**2.0 + disp(3*(i-1)+3)**2.0) ) maxdispnode = j
1165  maxdisp = max(maxdisp,sqrt(disp(3*(i-1)+1)**2.0 + disp(3*(i-1)+2)**2.0 + disp(3*(i-1)+3)**2.0))
1166  ENDIF
1167  ENDDO
1168 
1169 
1170 !
1171 ! Output results to a file
1172 !
1173  IF(mod(n,ifreq).EQ.0)THEN
1174  CALL mpi_reduce(maxdisp, gmaxdisp, 1, mpi_double_precision, mpi_max, 0, mpi_comm_world, ierr)
1175  IF(myid==0) WRITE(*,3000) n, t, gmaxdisp
1176  DO i=1,3*lnumnp
1177  WRITE(30,1001) local2global((i-1)/3+1),mod(i-1,3)+1,disp(i),v(i),a(i),t
1178  ENDDO
1179  WRITE(30,'()')
1180  !counter = 0
1181  !DO i = 1, NumNp
1182  ! IF(NodeProc(i) == myid) THEN
1183  ! counter = counter + 1
1184  ! !IF(Local2Global(i) == 41)THEN
1185  ! !IF(Local2Global(i) == 88) THEN
1186  ! IF(Local2Global(i) == 718) THEN
1187  ! IF(t<per1)THEN
1188  ! WRITE(40,*) t, disp(3*counter-2),disp(3*counter-1),disp(3*counter), SIN(3.14159*t/per1)
1189  ! ELSE
1190  ! WRITE(40,*) t, disp(3*counter-2),disp(3*counter-1),disp(3*counter), 0.0
1191  ! ENDIF
1192  ! EXIT
1193  ! ENDIF
1194  ! ENDIF
1195  !ENDDO
1196  ENDIF
1197 
1198  ENDDO
1199 
1200  ! Stop timing of program
1201  elapsedtime = elapsedtime + mpi_wtime()
1202  if(myid==0) print*,'Elapsed time of program: ',elapsedtime,' seconds'
1203 
1204  ! Finalize BS95
1205  CALL bs95free(bs95debug)
1206  CALL bs95finalize(bs95debug)
1207 
1208  ! Close output files
1209  CLOSE(30)
1210  !CLOSE(40)
1211 
1212  ! Deallocate arrays
1213  DEALLOCATE(cval_m)
1214  DEALLOCATE(aval_m)
1215  DEALLOCATE(rp_m)
1216  DEALLOCATE(cval_k)
1217  DEALLOCATE(aval_k)
1218  DEALLOCATE(rp_k)
1219  DEALLOCATE(cval_meff)
1220  DEALLOCATE(aval_meff)
1221  DEALLOCATE(rp_meff)
1222  DEALLOCATE(fext)
1223  DEALLOCATE(pbar)
1224  DEALLOCATE(fint)
1225  DEALLOCATE(disp)
1226  DEALLOCATE(v)
1227  DEALLOCATE(a)
1228  DEALLOCATE(s)
1229  DEALLOCATE(c)
1230  DEALLOCATE(ci)
1231  DEALLOCATE(local2global)
1232  DEALLOCATE(global2local)
1233  DEALLOCATE(nodeproc)
1234  IF (nprocs > 1) THEN
1235  DEALLOCATE(commprocs1)
1236  DEALLOCATE(numcommnodes1)
1237  DEALLOCATE(commnodes1)
1238  DEALLOCATE(commprocsfrom1)
1239  DEALLOCATE(numcommnodesfrom1)
1240  DEALLOCATE(commprocs2)
1241  DEALLOCATE(numcommnodes2)
1242  DEALLOCATE(commnodes2)
1243  DEALLOCATE(commprocsfrom2)
1244  DEALLOCATE(numcommnodesfrom2)
1245  ENDIF
1246 
1247  ! Formats
1248  1001 FORMAT(i5,' ',i1,' ',4(f15.6,' '))
1249  2000 FORMAT(100(f6.2,' '))
1250  3000 FORMAT(' ',i5,' ',f12.3,' ',f10.6)
1251 
1252 END PROGRAM fractography3d
1253 
initcomm1
subroutine initcomm1(global)
Definition: InitComm.f90:74
comp_row_global
Definition: comp_row_global.f90:26
removebcs_newa
subroutine removebcs_newa(nstart, ndim, a, newndim, newa)
Definition: removeBCs.f90:373
vol_elem_mat
subroutine vol_elem_mat(e, xnu, ci, cj, numat_vol, Integration)
Definition: vol_elem_mat.f90:53
initcomm2
subroutine initcomm2(global)
Definition: InitComm.f90:265
max
Vector_n max(const Array_n_const &v1, const Array_n_const &v2)
Definition: Vector_n.h:354
comp_row_add
subroutine comp_row_add(ndim, gndim, nrows1, nrows2, nnz1, nnz2, nstart1, rp1, cval1, aval1, nstart2, rp2, cval2, aval2)
Definition: comp_row_support.f90:89
Attribute::status
int status() const
Obtain the status of the attribute.
Definition: Attribute.h:240
removebcs_meff
subroutine removebcs_meff(ndim, nrows, nnz, nstart, rp1, cval, aval, newnrows, newnstart, newndim, global)
Definition: removeBCs.f90:84
sqrt
double sqrt(double d)
Definition: double.h:73
implicit_global
Definition: implicit_global.f90:76
feminp
subroutine feminp(glb, myid)
Definition: feminp.f90:53
Line_2::c
RT c() const
Definition: Line_2.h:150
comp_row_getval
subroutine comp_row_getval(ndim, nrows, nnz, nstart, rp, cval, aval, ipos, jpos, val)
Definition: comp_row_support.f90:511
fractography3d
program fractography3d
Definition: fractography3d.f90:101
implicit_bc_enforce
subroutine implicit_bc_enforce(NumNp, LocNumNp, disp, v, a, node_flag, boundary_value, t, myid)
Definition: implicit_bc_enforce.f90:83
comp_row_resize
subroutine comp_row_resize(ndim, nrows, nnz, nstart, rp, cval, aval, newnstart, newnrows)
Definition: comp_row_support.f90:592
implicit_v3d8_mass
subroutine implicit_v3d8_mass(NumEl, NumNP, NumMat, coor, nodes, MatType, ri, rho, xm, iElStart, iElEnd)
Definition: implicit_v3d8_mass.f90:82
get_fint
subroutine get_fint(ndim, nrows, nnz, nstart, rp_k, cval_k, aval_k, dispin, fint)
Definition: get_fint.f90:82
form
**********************************************************************Rocstar Simulation Suite Illinois Rocstar LLC All rights reserved ****Illinois Rocstar LLC IL **www illinoisrocstar com **sales illinoisrocstar com WITHOUT WARRANTY OF ANY **EXPRESS OR INCLUDING BUT NOT LIMITED TO THE WARRANTIES **OF FITNESS FOR A PARTICULAR PURPOSE AND **NONINFRINGEMENT IN NO EVENT SHALL THE CONTRIBUTORS OR **COPYRIGHT HOLDERS BE LIABLE FOR ANY DAMAGES OR OTHER WHETHER IN AN ACTION OF TORT OR **Arising OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE **USE OR OTHER DEALINGS WITH THE SOFTWARE **********************************************************************INTERFACE SUBROUTINE form
Definition: RFLO_ReadDataFileInt.h:26
removebcs_pbar
subroutine removebcs_pbar(nstart, ndim, pbar, newndim, newpbar)
Definition: removeBCs.f90:304
Aff_transformation_repS2::print
virtual std::ostream & print(std::ostream &os) const
Definition: Aff_transformationS2.h:277
min
Vector_n min(const Array_n_const &v1, const Array_n_const &v2)
Definition: Vector_n.h:346
implicit_v3d8_mass_consistent
subroutine implicit_v3d8_mass_consistent(NumEl, NumNP, NumMat, coor, nodes, MatType, ri, rho, ElConnVol)
Definition: implicit_v3d8_mass_consistent.f90:83
get_mat_stiffness
subroutine get_mat_stiffness(e, dnu, dmat)
Definition: get_mat_stiffness.f90:54
implicit_v3d8_me_k
subroutine implicit_v3d8_me_k(coor, ElConnVol, ci, numnp, nstart, nend, NumEL, MatType, NumMatType, enhanced_map, mixed_map, NumMatVol)
Definition: implicit_v3d8_me_K.f90:84
int_buf
Definition: fractography3d.f90:112
comp_row_addval
subroutine comp_row_addval(ndim, nrows, nnz, nstart, rp, cval, aval, ipos, jpos, val)
Definition: comp_row_support.f90:410