OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
nlsqf.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23C IDEBUG=1 -> output debug information into .out file
24#define IDEBUG 0
25C IOCSV=1 -> output fitting to curv.csv
26#define IOCSV 0
27!> @param[in] IFUNCS SUBROUTINE handle/ID for evaluation : 1=OGDEN0/OGDEN1
28!> @param[out] A updated coefficients
29!> @param[out] IRET return code
30!> 0 : optimization succeeded;
31!> 1 : optimization failed;
32!||====================================================================
33!|| mrqmin ../starter/source/materials/tools/nlsqf.F
34!||--- called by ------------------------------------------------------
35!|| law69_nlsqf ../starter/source/materials/tools/nlsqf.F
36!|| law69_nlsqf_auto ../starter/source/materials/tools/law69_nlsqf_auto.F
37!||--- calls -----------------------------------------------------
38!|| covsrt ../starter/source/materials/tools/nlsqf.F
39!|| inversion ../starter/source/materials/tools/nlsqf.F
40!|| mrqcof ../starter/source/materials/tools/nlsqf.F
41!||====================================================================
42 SUBROUTINE mrqmin(X,Y,SIG,NDATA,A,IA,
43 . MA,COVAR,ALPHA,NCA, ERRNOW,
44 . IFUNCS,GAMMA,IRET,ICOMP,MMAX,ATRY )
45C-----------------------------------------------
46C I m p l i c i t T y p e s
47C-----------------------------------------------
48#include "implicit_f.inc"
49C-----------------------------------------------
50C C o m m o n B l o c s
51C-----------------------------------------------
52#include "units_c.inc"
53C-----------------------------------------------
54C D u m m y A r g u m e n t s
55C-----------------------------------------------
56 INTEGER , INTENT(IN) :: MMAX
57 my_real , INTENT(INOUT) :: ATRY(MMAX)
58C----------------------------------------------------------------
59C L O C A L V A R I B L E S
60C----------------------------------------------------------------
61 INTEGER MA,NCA,NDATA,IA(MA) ,ICOMP
62 my_real
63 . gamma,errnow,a(ma),alpha(nca,nca),
64 . covar(nca,nca),x(ndata),y(ndata),sig(ndata)
65 INTEGER J,K,L,MFIT ,NMAX
66 PARAMETER (NMAX = 20) ! NMAX = MMAX
67 my_real errpre,beta(nmax),da(nmax)
68 SAVE errpre,beta,da,mfit
69 INTEGER IFUNCS, IRET
70C=======================================================================
71
72 iret = 0
73
74 IF (gamma < zero) THEN
75 mfit=0
76 DO j=1,ma
77 IF (ia(j)/=0) mfit=mfit+1
78 ENDDO
79! in Numerical Recript, 0.001 is used, however from my tests, the
80! difference is small (0.01 vs 0.001)
81 gamma=0.01
82 CALL mrqcof(x,y,sig,ndata,
83 . a,ia,ma,alpha,beta,nca,errnow,
84 . ifuncs,icomp)
85
86 errpre=errnow
87 DO j=1,ma
88 atry(j)=a(j)
89 ENDDO
90 ENDIF
91C
92 DO j=1,mfit
93 DO k=1,mfit
94 covar(j,k)=alpha(j,k)
95 ENDDO
96 covar(j,j)=alpha(j,j)*(1.+gamma)
97 da(j)=beta(j)
98 ENDDO
99C
100 CALL inversion(covar,mfit,nca,da,1,1,iret)
101 IF (iret /= 0) THEN
102 IF (idebug > 0) THEN
103 WRITE(*,*) ' IRET = ', iret
104 DO j=1,ma
105 WRITE(iout, *) 'J, A(J) = ', j, a(j)
106 ENDDO
107 ENDIF
108 RETURN
109 ENDIF
110
111 IF (gamma == zero) THEN
112 CALL covsrt(covar,nca,ma,ia,mfit)
113 RETURN
114 ENDIF
115C
116 j=0
117 DO l=1,ma
118 IF(ia(l)/=0) THEN
119 j=j+1
120 atry(l)=a(l)+da(j)
121 ENDIF
122 ENDDO
123
124 IF (idebug > 0) THEN
125 WRITE(iout,*) __file__,__line__
126 DO j=1, ma
127 write(iout,*) 'J,ATRY(J) = ', j, atry(j)
128 ENDDO
129 ENDIF
130
131c limit ATRY (ALPHA) to avoid too big power (overflow)
132c |ALPHA_i| <= 20.0
133 DO j=1, ma-1, 2
134 if (atry(j+1) > twenty) then
135 IF (idebug > 0) THEN
136 write(iout,*) 'Setting Big ALPHA to ',twenty
137 write(iout,*) ' J, ATRY(J) = ',j+1,atry(j+1)
138 ENDIF
139 atry(j+1) = twenty
140 ELSEIF (atry(j+1) < -twenty) THEN
141 IF (idebug > 0) THEN
142 write(iout,*) 'Setting Big ALPHA to -',twenty
143 write(iout,*) ' J, ATRY(J) = ',j+1,atry(j+1)
144 ENDIF
145 atry(j+1) = -twenty
146 ENDIF
147 ENDDO
148
149C
150 CALL mrqcof(x,y,sig,ndata,
151 . atry,ia,ma,covar,da,nca,errnow,
152 . ifuncs,icomp)
153C
154 IF (errnow < errpre) THEN
155 gamma=0.1*gamma
156 errpre=errnow
157 DO j=1,mfit
158 DO k=1,mfit
159 alpha(j,k)=covar(j,k)
160 ENDDO
161 beta(j)=da(j)
162 ENDDO
163 DO l=1,ma
164 a(l)=atry(l)
165 ENDDO
166 ELSE
167 gamma=10.*gamma
168 errnow=errpre
169 ENDIF
170C-----------
171 RETURN
172 END SUBROUTINE mrqmin
173C-----------
174C
175C
176!||====================================================================
177!|| mrqcof ../starter/source/materials/tools/nlsqf.F
178!||--- called by ------------------------------------------------------
179!|| mrqmin ../starter/source/materials/tools/nlsqf.F
180!||--- calls -----------------------------------------------------
181!|| my_exit ../starter/source/output/analyse/analyse.c
182!|| ogden0 ../starter/source/materials/tools/nlsqf.F
183!|| ogden1 ../starter/source/materials/tools/nlsqf.F
184!||====================================================================
185 SUBROUTINE mrqcof(X,Y,SIG, NDATA,
186 * A,IA,MA,ALPHA,BETA,NALP,
187 * ERRNOW,IFUNCS,ICOMP)
188C-----------------------------------------------
189C I m p l i c i t T y p e s
190C-----------------------------------------------
191#include "implicit_f.inc"
192C-----------------------------------------------
193C D u m m y A r g u m e n t s
194C-----------------------------------------------
195 INTEGER MA,NALP,NDATA,IA(MA),ICOMP
196 my_real
197 . ERRNOW,A(MA),ALPHA(NALP,NALP),BETA(MA),
198 . X(NDATA),Y(NDATA),SIG(NDATA)
199 INTEGER IFUNCS
200C----------------------------------------------------------------
201C L O C A L V A R I B L E S
202C----------------------------------------------------------------
203 INTEGER MFIT,I,J,K,L,M,MMAX
204 PARAMETER (MMAX=20)
205 my_real
206 . dy,wt,ymod,dyda(mmax)
207
208 my_real y_min
209
210 y_min = em3
211
212C=======================================================================
213 mfit=0
214 DO 11 j=1,ma
215 IF (ia(j)/=0) mfit=mfit+1
21611 CONTINUE
217 DO 13 j=1,mfit
218 DO 12 k=1,j
219 alpha(j,k)=0.
22012 CONTINUE
221 beta(j)=0.
22213 CONTINUE
223 errnow=zero
224 IF(icomp == 0) THEN
225 DO 16 i=1,ndata
226 IF (ifuncs == 1) THEN
227 CALL ogden0(x(i),a,ymod, ma)
228 CALL ogden1(x(i),a,dyda, ia, ma)
229 ELSE
230! TBD, UNKOWN IFUNCS, ERROR OUT!!!
231 WRITE(*,*) 'ERROR: UNKOWN IFUNCS = ', ifuncs
232 CALL my_exit(2)
233 ENDIF
234
235 dy=y(i)-ymod
236 errnow=errnow+dy*dy/(sig(i)*sig(i))
237c IF (Y(I)>EM5) DY=Y(I)-YMOD
238c IF (Y(I)<EM5) DY=EM5
239
240 j=0
241 DO 15 l=1,ma
242 IF(ia(l)/=0) THEN
243 j=j+1
244 wt=dyda(l)
245 wt=dyda(l)/(sig(i)*sig(i))
246 k=0
247 DO 14 m=1,l
248 IF(ia(m)/=0) THEN
249 k=k+1
250 alpha(j,k)=alpha(j,k)+wt*dyda(m)
251 ENDIF
25214 CONTINUE
253 beta(j)=beta(j)+dy*wt
254 ENDIF
25515 CONTINUE
25616 CONTINUE
257 ELSE
258 DO i=1,ndata
259 IF (ifuncs == 1) THEN
260 CALL ogden0(x(i),a,ymod, ma)
261 CALL ogden1(x(i),a,dyda, ia, ma)
262 ELSE
263! TBD, UNKOWN IFUNCS, ERROR OUT!!!
264 WRITE(*,*) 'ERROR: UNKOWN IFUNCS = ', ifuncs
265 CALL my_exit(2)
266 ENDIF
267
268 dy=y(i)-ymod
269 errnow = errnow + dy*dy
270 j=0
271 DO l=1,ma
272 IF(ia(l)/=0) THEN
273 j=j+1
274 wt=dyda(l)
275 k=0
276 DO m=1,l
277 IF(ia(m)/=0) THEN
278 k=k+1
279 alpha(j,k)=alpha(j,k)+wt*dyda(m)
280 ENDIF
281 ENDDO ! M
282 beta(j)=beta(j)+dy*wt
283 ENDIF
284 ENDDO ! L
285 ENDDO ! I
286 ENDIF
287 DO 18 j=2,mfit
288 DO 17 k=1,j-1
289 alpha(k,j)=alpha(j,k)
29017 CONTINUE
291
29218 CONTINUE
293 RETURN
294 END SUBROUTINE mrqcof
295C
296
297!||====================================================================
298!|| inversion ../starter/source/materials/tools/nlsqf.F
299!||--- called by ------------------------------------------------------
300!|| mrqmin ../starter/source/materials/tools/nlsqf.F
301!|| mrqmin_law92 ../starter/source/materials/mat/mat092/law92_nlsqf.F90
302!||====================================================================
303 SUBROUTINE inversion(A,N,NP,B,M,MP,IRET)
304C-----------------------------------------------
305C I M P L I C I T T Y P E S
306C-----------------------------------------------
307#include "implicit_f.inc"
308C-----------------------------------------------
309C D u m m y A r g u m e n t s
310C-----------------------------------------------
311 INTEGER M,MP,N,NP
312 my_real A(NP,NP),B(NP,MP)
313C----------------------------------------------------------------
314C L O C A L V A R I B L E S
315C----------------------------------------------------------------
316 INTEGER I,NMAX,ICOL,IROW,J,K,L,LL
317 PARAMETER (NMAX=50)
318 integer indxc(nmax),indxr(nmax),ipiv(nmax)
319 my_real
320 . big,dum,pivinv
321
322 INTEGER IRET
323
324 iret = 0
325
326 icol = -1
327 irow = -1
328
329C=======================================================================
330 DO 11 j=1,n
331 ipiv(j)=0
33211 CONTINUE
333 DO 22 i=1,n
334 big=zero
335 DO 13 j=1,n
336 IF(ipiv(j)/=1)THEN
337 DO 12 k=1,n
338 IF (ipiv(k)==0) THEN
339 IF (abs(a(j,k))>=big)THEN
340 big=abs(a(j,k))
341 irow=j
342 icol=k
343 ENDIF
344 ELSE IF (ipiv(k)>1) THEN
345 iret = __line__
346 WRITE(*,*) 'SINGULAR MATRIX IN INVERSION'
347 RETURN
348 ENDIF
34912 CONTINUE
350 ENDIF
35113 CONTINUE
352
353 IF (icol == -1) THEN
354 iret = __line__
355 WRITE(*,*) 'FATAL ERROR IN INVERSION'
356 RETURN
357 ENDIF
358
359 ipiv(icol)=ipiv(icol)+1
360 IF (irow/=icol) THEN
361 DO 14 l=1,n
362 dum=a(irow,l)
363 a(irow,l)=a(icol,l)
364 a(icol,l)=dum
36514 CONTINUE
366 DO 15 l=1,m
367 dum=b(irow,l)
368 b(irow,l)=b(icol,l)
369 b(icol,l)=dum
37015 CONTINUE
371 ENDIF
372 indxr(i)=irow
373 indxc(i)=icol
374 IF (a(icol,icol)==0.) THEN
375 iret = __line__
376 WRITE(*,*) 'SINGULAR MATRIX IN INVERSION'
377 RETURN
378 ENDIF
379 pivinv=1./a(icol,icol)
380 a(icol,icol)=1.
381 DO 16 l=1,n
382 a(icol,l)=a(icol,l)*pivinv
38316 CONTINUE
384 DO 17 l=1,m
385 b(icol,l)=b(icol,l)*pivinv
38617 CONTINUE
387 DO 21 ll=1,n
388 IF(ll/=icol)THEN
389 dum=a(ll,icol)
390 a(ll,icol)=0.
391 DO 18 l=1,n
392 a(ll,l)=a(ll,l)-a(icol,l)*dum
39318 CONTINUE
394 DO 19 l=1,m
395 b(ll,l)=b(ll,l)-b(icol,l)*dum
39619 CONTINUE
397 ENDIF
39821 CONTINUE
39922 CONTINUE
400 DO 24 l=n,1,-1
401 IF(indxr(l)/=indxc(l))THEN
402 DO 23 k=1,n
403 dum=a(k,indxr(l))
404 a(k,indxr(l))=a(k,indxc(l))
405 a(k,indxc(l))=dum
40623 CONTINUE
407 ENDIF
40824 CONTINUE
409C-----------
410 RETURN
411 END SUBROUTINE inversion
412C
413C
414!||====================================================================
415!|| covsrt ../starter/source/materials/tools/nlsqf.F
416!||--- called by ------------------------------------------------------
417!|| mrqmin ../starter/source/materials/tools/nlsqf.F
418!||====================================================================
419 SUBROUTINE covsrt(COVAR,SIZCOVAR,MA,IA,MFIT)
420C-----------------------------------------------
421C I M P L I C I T T Y P E S
422C-----------------------------------------------
423#include "implicit_f.inc"
424C----------------------------------------------------------------
425C D u m m y A r g u m e n t s
426C----------------------------------------------------------------
427 INTEGER MA,MFIT,SIZCOVAR,IA(MA)
428 my_real COVAR(SIZCOVAR,SIZCOVAR),SWAP
429C----------------------------------------------------------------
430C L O C A L V A R I B L E S
431C----------------------------------------------------------------
432 INTEGER I,J,K
433C=======================================================================
434 DO 12 I=mfit+1,ma
435 DO 11 j=1,i
436 covar(i,j)=0.
437 covar(j,i)=0.
43811 CONTINUE
43912 CONTINUE
440 k=mfit
441 DO 15 j=ma,1,-1
442 IF(ia(j)/=0)THEN
443 DO 13 i=1,ma
444 swap=covar(i,k)
445 covar(i,k)=covar(i,j)
446 covar(i,j)=swap
44713 CONTINUE
448 DO 14 i=1,ma
449 swap=covar(k,i)
450 covar(k,i)=covar(j,i)
451 covar(j,i)=swap
45214 CONTINUE
453 k=k-1
454 ENDIF
45515 CONTINUE
456C-----------
457 RETURN
458 END SUBROUTINE covsrt
459
460!> @purpose
461!> calculate stress from stretch for OGDEN (uniaxial tension)
462!> @param[in] X stretch
463!> @param[in] A(NA) material parameters [mu_1, alpha_1,mu_2, alpha_2,...]
464!> @param[out] Y stress
465!||====================================================================
466!|| ogden0 ../starter/source/materials/tools/nlsqf.F
467!||--- called by ------------------------------------------------------
468!|| law69_nlsqf ../starter/source/materials/tools/nlsqf.F
469!|| law69_nlsqf_auto ../starter/source/materials/tools/law69_nlsqf_auto.F
470!|| mrqcof ../starter/source/materials/tools/nlsqf.F
471!||====================================================================
472 SUBROUTINE ogden0(X,A,Y,NA)
473C-----------------------------------------------
474C I m p l i c i t T y p e s
475C-----------------------------------------------
476#include "implicit_f.inc"
477C-----------------------------------------------
478C D u m m y A r g u m e n t s
479C-----------------------------------------------
480 INTEGER NA
481 my_real X, A(NA), Y
482C----------------------------------------------------------------
483C L O C A L V A R I B L E S
484C----------------------------------------------------------------
485 INTEGER I
486C=======================================================================
487 Y=zero
488 DO i=1,na-1,2
489 y = y + a(i)*(x**(-1.+a(i+1))-x**(-0.5*a(i+1)-1.))
490 ENDDO
491C-----------
492 RETURN
493 END SUBROUTINE ogden0
494
495!> @purpose
496!> calculate d(stress)/d(A) OGDEN (uniaxial tension)
497!> @param[in] X stretch
498!> @param[in] A(NA) material parameters [mu_1, alpha_1,mu_2, alpha_2,...]
499!> @param[in] IA(NA) flags if derivative needs to be calculated
500!> 1=yes, 0=no
501!> @param[out] DYDA(NA) derivatives d(stress)/d(A)
502!||====================================================================
503!|| ogden1 ../starter/source/materials/tools/nlsqf.F
504!||--- called by ------------------------------------------------------
505!|| mrqcof ../starter/source/materials/tools/nlsqf.F
506!||====================================================================
507 SUBROUTINE ogden1(X,A,DYDA,IA,NA)
508C-----------------------------------------------
509C I m p l i c i t T y p e s
510C-----------------------------------------------
511#include "implicit_f.inc"
512C-----------------------------------------------
513C D u m m y A r g u m e n t s
514C-----------------------------------------------
515 INTEGER NA
516 INTEGER IA(NA)
517 my_real X, A(NA),DYDA(NA)
518 my_real TMP1, TMP2
519C----------------------------------------------------------------
520C L O C A L V A R I B L E S
521C----------------------------------------------------------------
522 INTEGER I
523C=======================================================================
524 DO I=1,na-1,2
525 tmp1 = x**(-1.+a(i+1))
526 tmp2 = x**(-0.5*a(i+1)-1.)
527 IF (ia(i) /= 0) THEN
528 dyda(i) = tmp1 - tmp2
529 ENDIF
530
531 IF (ia(i+1) /= 0) THEN
532 dyda(i+1) = a(i)*log(x)*(tmp1 + 0.5*tmp2)
533 ENDIF
534 ENDDO
535
536 RETURN
537 END SUBROUTINE ogden1
538
539!> @param[in] LAWID law_ID on /MAT/LAW69 (1=Ogden, 2=Mooney-Rivlin)
540!> @param[in] ICHECK checking level
541!> <=0 no constraint (IVALID always return as 1)
542!> 1 SUM( mu(i) * alpha(i) ) > 0.0
543!> 2 mu(i) * alpha(i) > 0.0
544!> 3 Try ICHECK=2 at first, if it fails in fitting, switch to ICHECK=1, and try again
545!> @param[in] ERRTOL If ERRAVE < ERRTOL, data fitting converges.
546!> ERRAVE = ( SUM [ ABS ( ( Y_inp-Y_fit) / Y_inp ) ) / NPT
547!> @param[out] MUAL(1:NMUAL) optimized material properties
548!||====================================================================
549!|| law69_nlsqf ../starter/source/materials/tools/nlsqf.F
550!||--- called by ------------------------------------------------------
551!|| law69_upd ../starter/source/materials/mat/mat069/law69_upd.F
552!||--- calls -----------------------------------------------------
553!|| ancmsg ../starter/source/output/message/message.F
554!|| law69_check ../starter/source/materials/tools/nlsqf.F
555!|| law69_guess1 ../starter/source/materials/tools/nlsqf.F
556!|| law69_guess_bounds ../starter/source/materials/tools/nlsqf.F
557!|| mrqmin ../starter/source/materials/tools/nlsqf.F
558!|| my_exit ../starter/source/output/analyse/analyse.c
559!|| ogden0 ../starter/source/materials/tools/nlsqf.F
560!||--- uses -----------------------------------------------------
561!|| message_mod ../starter/share/message_module/message_mod.F
562!||====================================================================
563 SUBROUTINE law69_nlsqf(LAWID,STRETCH,Y,MA,NMUAL,NPT,MUAL,
564 $ ICHECK, NSTART, ERRTOL,ID,TITR,IS_ENCRYPTED)
565 USE message_mod
566 USE names_and_titles_mod , ONLY : nchartitle
567C-----------------------------------------------
568C I m p l i c i t T y p e s
569C-----------------------------------------------
570#include "implicit_f.inc"
571C-----------------------------------------------
572C C o m m o n B l o c k s
573C-----------------------------------------------
574#include "units_c.inc"
575C-----------------------------------------------
576 INTEGER MAXA
577 PARAMETER (MAXA=20)
578
579 integer lawid, npt, ma, icheck, nstart
580 INTEGER , INTENT(IN) :: IS_ENCRYPTED
581 my_real errtol
582 INTEGER I,NONZERO(MAXA),ITER,J,NMUAL
583 my_real GAMMA,ERRNOW,ERRPRE,
584 . stretch(npt),mual(10)
585 my_real a(maxa),covar(maxa,maxa),alpha(maxa,maxa),y(npt)
586 my_real mcof_min(maxa), mcof_max(maxa)
587 my_real yogd
588
589 INTEGER ID
590 CHARACTER(LEN=NCHARTITLE) :: TITR
591 INTEGER MINITER_LM, MAXITER_LM
592 DATA miniter_lm /5/
593 SAVE miniter_lm
594
595 DATA maxiter_lm /100/
596 SAVE maxiter_lm
597
598 ! if ( (ERRNOW-ERRPRE) / ERRPRE ) < EPS_LM, CNT_HIT_EPS_LM = CNT_HIT_EPS_LM + 1,
599 ! if CNT_HIT_EPS_LM == LMT_HIT_EPS_LM, it converges
600 my_real eps_lm
601 DATA eps_lm /1e-3/
602 SAVE eps_lm
603
604 INTEGER CNT_HIT_EPS_LM
605 INTEGER LMT_HIT_EPS_LM
606 DATA lmt_hit_eps_lm /2/
607 SAVE lmt_hit_eps_lm
608
609 INTEGER LMSTOP
610
611 INTEGER IFUNCS
612
613 INTEGER NGUESS
614 my_real a0(maxa)
615 my_real errave, errave_min, err
616
617 INTEGER ISTART, NPSAVED, IVALID,ICOMP, MMAX
618 parameter(mmax =20)
619 my_real atry(mmax)
620
621 INTEGER ICURV
622 LOGICAL lopened
623
624c during LM optimization, if the objective is not improved,
625c GAMMA will be increased. We should terminate the iteration once
626c GAMMA becomes very huge.
627 my_real gamma_stop
628 data gamma_stop /1e15/
629 save gamma_stop
630
631c if check the validity of initial guess (0=no, 1=yes)
632c we don't want to check, because invalid initial point could converge to valid point.
633 INTEGER ICHECK_GUESS
634 data icheck_guess /0/
635 save icheck_guess
636
637 INTEGER JCHECK, IFIT_SUCCESS
638
639c if enforce mu(i) < mu(i+1) in generating initial guess
640c we don't need this anymore since we are using random numbers instead of loop through all
641c the combinations
642 INTEGER MU_INCR_GUESS
643 data MU_INCR_GUESS /0/
644 save MU_INCR_GUESS
645
646 my_real MAX_ABS_YI
647 my_real small_fac_abs_yi, small_abs_yi
648
649 my_real, DIMENSION(:), ALLOCATABLE :: sig
650 my_real spready
651 INTEGER IRET
652 INTEGER IRND1
653! ------------------------------------------------------------------------------
654 ALLOCATE (sig(1:npt))
655 errave_min = zero
656!
657 !
658 IF ( maxa < ma ) THEN
659 WRITE(*,*) 'ERROR, MAXA < MA'
660 WRITE(*,*) __file__,__line__
661 CALL my_exit(2)
662 ENDIF
663 icomp=0
664 IF(icheck < 0) icomp= 1
665 icheck = abs(icheck)
666 ! IF ABS(Y(I)) < SMALL_ABS_YI, use SMALL_ABS_YI to avoid
667 ! unnecessary numerical issues when avoid divided by small value
668
669 max_abs_yi = zero
670 DO i = 1, npt
671 max_abs_yi = max( max_abs_yi, abs(y(i)) )
672 ENDDO
673
674 small_fac_abs_yi = em3
675
676 small_abs_yi = max_abs_yi * small_fac_abs_yi
677
678 IF (idebug > 0) THEN
679 WRITE(iout, *) ' MAX_ABS_YI, SMALL_FAC_ABS_YI, SMALL_ABS_YI = ',
680 $ max_abs_yi, small_fac_abs_yi, small_abs_yi
681 ENDIF
682
683 spready=0.01
684 IF(icomp == 0) THEN
685 DO j=1,npt
686 sig(j)=spready*max(small_abs_yi, abs(y(j)) )
687 IF (idebug > 0) THEN
688 WRITE(iout, *) 'J, SIG(J) = ', j, sig(j)
689 ENDIF
690 ENDDO
691 ELSE
692 DO j=1,npt
693 sig(j) = y(j)
694 IF(sig(j) == zero) sig(j) = em15
695 IF (idebug > 0) THEN
696 WRITE(iout, *) 'J, SIG(J) = ', j, sig(j)
697 ENDIF
698 ENDDO
699 ENDIF
700
701 IF (maxa < ma) THEN
702! TBD_KANG
703 WRITE(*,*) 'ERROR: MAXA < MA'
704 WRITE(*,*) ' MAXA = ', maxa
705 WRITE(*,*) ' MA = ', maxa
706 WRITE(*,*) ' FILE = ', __file__
707 WRITE(*,*) ' LINE = ', __line__
708 CALL my_exit(2)
709 ENDIF
710
711C=======================================================================
712 DO i=1,nmual
713 nonzero(i)=1
714 ENDDO
715
716 IF (lawid == 1) then ! OGDEN
717 DO i=nmual+1,10
718 nonzero(i)=0
719 a(i)=zero
720 ENDDO
721
722 ifuncs = 1 ! SUBROUTINE OGDEN0/OGDEN1
723
724 ELSEIF (lawid == 2) then ! Mooney-Rivlin
725 DO i=nmual+1,ma
726 nonzero(i)=0
727 ENDDO
728 nonzero(2)=0
729 nonzero(4)=0
730
731 DO i=nmual+1,ma
732 a(i)=zero
733 ENDDO
734
735 ifuncs = 1 ! SUBROUTINE OGDEN0/OGDEN1
736
737 ENDIF
738
739 jcheck =icheck
740 IF (jcheck == 3) THEN
741 jcheck = 2
742 ENDIF
743
74499 CONTINUE ! location to start optimization with different JCHECK
745
746 ifit_success = 0
747
748 istart = 0
749 npsaved = 0
750
751c calculate MCOF_MIN, MCOF_MAX
752 CALL law69_guess_bounds(lawid, jcheck, stretch, y, npt,
753 $ nmual, mcof_min, mcof_max,icomp)
754
755c initialize random number generator in LAW69_GUESS1
756 irnd1 = 205
757 atry(1:mmax) = zero ! Initialization of ATRY
758 DO 111 WHILE ( istart < nstart )
759c Get One Guess Point A0 (1: Nmual)
760 CALL law69_guess1(
761 $ lawid, nmual, icheck_guess, mu_incr_guess,irnd1,
762 $ a0, nonzero, mcof_min, mcof_max, nguess)
763
764 IF (nguess == 0) THEN
765 GOTO 112
766 ENDIF
767
768C calculate averaged ERROR before LM
769 errave = 0.0
770 IF(icomp == 0) THEN
771 DO i=1,npt
772 CALL ogden0(stretch(i), a0, yogd, nmual)
773 err = abs(y(i) - yogd) / max(small_abs_yi, abs(y(i)))
774 errave = errave + err
775 ENDDO
776 ELSE
777 DO i=1,npt
778 CALL ogden0(stretch(i), a0, yogd, nmual)
779 err = abs(y(i) - yogd) / max(em15,abs(y(i)))
780 errave = errave + err
781 ENDDO
782 ENDIF
783C
784 errave = errave / (1.0 * npt)
785
786 istart = istart + 1
787
788 IF (idebug > 0) THEN
789 WRITE(iout,*) ' ISTART = ', istart
790 WRITE(iout,*) ' Before LM optimization ...'
791 DO i = 1, nmual
792 WRITE(iout, *) ' I, A0(I) = ', i, a0(i)
793 ENDDO
794 WRITE(iout,*) ' LM0, ERRAVE = ', errave
795 ENDIF
796
797! start loop of LM optimization
798 iter = 0
799 cnt_hit_eps_lm = 0
800
801 iter = iter + 1
802 gamma=-1
803 CALL mrqmin(stretch,y,sig,npt,a0,nonzero,
804 $ nmual,covar,alpha,ma,errnow,
805 $ ifuncs,gamma,iret,icomp,mmax,atry)
806 IF (iret /= 0) GOTO 111
807
808 IF (idebug > 0) THEN
809 errave = 0.0
810 IF(icomp == 0) THEN
811 DO i=1,npt
812 CALL ogden0(stretch(i), a0, yogd, nmual)
813 err = abs(y(i) - yogd) / max(small_abs_yi, abs(y(i)))
814 errave = errave + err
815 ENDDO
816 ELSE
817 DO i=1,npt
818 CALL ogden0(stretch(i), a0, yogd, nmual)
819 err = abs(y(i) - yogd) / max(em15,abs(y(i)))
820 errave = errave + err
821 ENDDO
822 ENDIF
823 errave = errave / (1.0 * npt)
824
825 WRITE(iout, '(A,I4, 3E16.8)')
826 $ 'ITER, ERRNOW, GAMMA, ERRAVE = ',
827 $ iter, errnow, gamma, errave
828 DO i = 1, nmual
829 WRITE(iout, *) ' - I, A0(I) = ', i, a0(i)
830 ENDDO
831 ENDIF
832
83321 CONTINUE
834 errpre=errnow
835
836 iter = iter + 1
837 CALL mrqmin(stretch,y,sig,npt,a0,nonzero,
838 $ ma,covar,alpha, ma, errnow,
839 $ ifuncs,gamma,iret,icomp,mmax,atry)
840 IF (iret /= 0) GOTO 111
841
842 IF (idebug == 1) THEN
843 errave = zero
844 IF(icomp == 0) THEN
845 DO i=1,npt
846 CALL ogden0(stretch(i), a0, yogd, nmual)
847 err = abs(y(i) - yogd) / max(small_abs_yi, abs(y(i)))
848 errave = errave + err
849 ENDDO
850 ELSE
851 DO i=1,npt
852 CALL ogden0(stretch(i), a0, yogd, nmual)
853 err = abs(y(i) - yogd) / max(em15,abs(y(i)))
854 errave = errave + err
855 ENDDO
856 ENDIF
857 errave = errave / (1.0 * npt)
858
859 WRITE(iout, '(A,I4, 5E16.8)')
860 $ 'ITER, ERRNOW, GAMMA, ERRAVE, ERRNOW-ERRPRE,'//
861 $ '(ERRNOW-ERRPRE)/ERRPRE = ',
862 $ iter, errnow, gamma, errave, errnow-errpre,
863 $ (errnow-errpre)/errpre
864 DO i = 1, nmual
865 WRITE(iout, *) ' - I, A0(I) = ', i, a0(i)
866 ENDDO
867 ENDIF
868
869c IF A0(J) is too small, restart from next initial guess
870 DO j = 1, nmual
871 IF ( abs( a0(j) ) < 1e-20 ) THEN
872 goto 111 ! restart from next initial guess
873c WRITE(IOUT, *) ' A is too small, and could result in error'
874c WRITE(IOUT,*) 'J, A0(J) = ', J, A0(J)
875c STOP
876 ENDIF
877 ENDDO
878
879c check convergence of LM optimization
880 lmstop = 0
881 IF (idebug > 0) THEN
882 WRITE(iout, *) ' ERRNOW/(1.0*NPT) = ', errnow/(1.0*npt)
883 WRITE(iout, *) ' ERRNOW < ERRPRE = ', (errnow < errpre)
884 ENDIF
885
886 IF ( iter > miniter_lm ) THEN
887 IF (errnow < errpre) THEN
888 IF ( abs(errpre) > zero) THEN
889 IF ( abs( (errnow-errpre)/ errpre ) < eps_lm) THEN
890 cnt_hit_eps_lm = cnt_hit_eps_lm + 1
891
892 IF (idebug > 0) THEN
893 WRITE(iout,*)
894 $ ' CNT_HIT_EPS_LM, ABS((ERRNOW-ERRPRE)/ERRPRE) = ',
895 $ cnt_hit_eps_lm, abs( (errnow-errpre)/ errpre )
896 ENDIF
897
898 IF ( cnt_hit_eps_lm >= lmt_hit_eps_lm ) THEN
899 IF (idebug > 0) THEN
900 WRITE(iout,*) 'STOP AT ', __line__
901 ENDIF
902 lmstop = 1
903 ENDIF
904 ENDIF
905 ENDIF
906 ELSEIF (iter >= maxiter_lm .OR. gamma >= gamma_stop) THEN
907 IF (idebug > 0) THEN
908 WRITE(iout,*) 'STOP AT ', __line__
909 ENDIF
910 lmstop = 1
911 ENDIF
912 ENDIF
913
914 IF (lmstop== 0) THEN
915 GOTO 21
916 ENDIF
917! end loop of LM optimization
918
919 errave = 0.0
920 IF(icomp == 0) THEN
921 DO i=1,npt
922 CALL ogden0(stretch(i), a0, yogd, nmual)
923 err = abs(y(i) - yogd) / max(small_abs_yi, abs(y(i)))
924 errave = errave + err
925 ENDDO
926 ELSE
927 DO i=1,npt
928 CALL ogden0(stretch(i), a0, yogd, nmual)
929 err = abs(y(i) - yogd) / max(em15,abs(y(i)))
930 errave = errave + err
931 ENDDO
932 ENDIF
933 errave = errave / (1.0 * npt)
934
935 IF (idebug > 0) THEN
936 WRITE(iout,*) ' After LM optimization ...'
937 DO i = 1, ma
938 WRITE(iout, *) ' I, A0(I) = ', i, a0(i)
939 ENDDO
940 WRITE(iout,*) ' LM1, ERRAVE = ', errave
941 ENDIF
942
943 CALL law69_check(lawid, a0, nmual, jcheck, 0, ivalid)
944 IF (ivalid > 0) THEN
945 IF ( npsaved==0 .OR.
946 $ ( npsaved>0 .AND. errave<errave_min) ) THEN
947 npsaved = npsaved + 1
948 errave_min = errave
949 DO i = 1, nmual
950 a(i) = a0(i)
951 ENDDO
952 ENDIF
953 ELSE
954 IF (idebug > 0) THEN
955 WRITE(iout,*) __file__,__line__
956 WRITE(iout,*) ' LM converges to invalid point'
957 ENDIF
958 ENDIF
959
960 IF (npsaved > 0) THEN
961 IF (idebug > 0) THEN
962 WRITE(*, *) ' ISTART, NPSAVED, ERRAVE_MIN = ',
963 $ istart, npsaved, errave_min
964 WRITE(iout, *) ' ISTART, NPSAVED, ERRAVE_MIN = ',
965 $ istart, npsaved, errave_min
966 ENDIF
967
968 IF ( errave_min < errtol ) THEN
969 ifit_success = 1
970 GOTO 112
971 ENDIF
972 ELSE
973 IF (idebug > 0) THEN
974 WRITE(*, *) ' ISTART, NPSAVED ',
975 $ istart, npsaved
976 WRITE(iout, *) ' ISTART, NPSAVED ',
977 $ istart, npsaved
978 ENDIF
979 ENDIF
980
981111 CONTINUE ! WHILE ( ISTART < NSTART )
982
983112 CONTINUE
984
985 IF (ifit_success == 0 .AND. icheck == 3) THEN
986 IF (jcheck == 2) THEN
987 jcheck = 1
988 IF (idebug > 0) THEN
989 IF (npsaved > 0) THEN
990 WRITE(*,'(A)')
991 $ '*** Curve fitting result of /MAT/LAW69 with'
992 $ //' ICHECK=2 is not satisfactory.'
993 WRITE(*,'(A,F10.2,A)') 'ERRAVE = ', errave_min*100.0,'%'
994 WRITE(*,'(A)') ' Switching to ICHECK=1 and trying again.'
995
996 ELSE
997 WRITE(*,'(A)')
998 $ '*** Failed to fit /MAT/LAW69 with ICHECK=2.'
999 WRITE(*,'(A)') ' Switching to ICHECK=1 and trying again.'
1000 ENDIF
1001 ENDIF
1002 GOTO 99
1003 ENDIF
1004 ENDIF
1005
1006 DEALLOCATE (sig)
1007
1008 IF (ifit_success == 0) THEN
1009 IF (npsaved == 0) THEN
1010 CALL ancmsg(msgid=901,
1011 . msgtype=msgerror,
1012 . anmode=aninfo,
1013 . i1=id,
1014 . c1=titr)
1015 ENDIF
1016 ENDIF
1017C
1018 DO i=1,10
1019 mual(i)=a(i)
1020 ENDDO
1021
1022 IF (lawid == 2) THEN
1023 DO i=5,10
1024 mual(i) = zero
1025 ENDDO
1026 ENDIF
1027 IF(is_encrypted == 0)THEN
1028 WRITE(iout,'(//6X,A,/)')'FITTING RESULT COMPARISON:'
1029 WRITE(iout,'(6X,A,/)')'UNIAXIAL TEST DATA'
1030 WRITE(iout,'(A20,5X,A20,A30)') 'NOMINAL STRAIN',
1031 * 'NOMINAL STRESS(TEST)', 'NOMINAL STRESS(RADIOSS)'
1032 ENDIF
1033
1034c output curcves to .csv format to simplify visualization
1035 icurv = 0
1036 IF (iocsv> 0) THEN
1037 DO icurv = 25, 99
1038 inquire (unit=icurv, opened=lopened)
1039 if (.not. lopened) goto 77
1040 ENDDO
104177 CONTINUE
1042 OPEN(unit=icurv, file='curv.csv')
1043 WRITE(icurv,'(A)') 'NOMINAL STRAIN, NOMINAL STRESS(TEST), '//
1044 $ 'NOMINAL STRESS(RADIOSS)'
1045 ENDIF
1046 IF(is_encrypted == 0)THEN
1047 DO i=1,npt
1048 CALL ogden0(stretch(i), a, yogd, nmual)
1049!! WRITE(IOUT,'(F18.4,F20.4,F28.4)')
1050 WRITE(iout, '(1F18.6, 3X,1F20.13, 6X, 1F20.13)')
1051 * stretch(i)-one,y(i),yogd
1052 IF (icurv > 0) THEN
1053 WRITE(icurv, '(F18.4, A, F18.4, A, F18.4)')
1054 $ stretch(i)-one,',',y(i),',',yogd
1055 ENDIF
1056 ENDDO
1057 ENDIF
1058
1059 WRITE(iout, *) ''
1060 WRITE(iout, '(A)') '-------------------------------------------'
1061 WRITE(iout, '(A,F10.2,A)') 'AVERAGED ERROR OF FITTING : ',
1062 $ errave_min*100.0, '%'
1063
1064c output curcves to .csv format to simplify visualization
1065 IF ( icurv > 0) THEN
1066 CLOSE(icurv)
1067 ENDIF
1068C-----------
1069 RETURN
1070 END
1071
1072c if abs(value) < small, it means value is very close to zero,
1073c reset it to be small or -small to avoid potential numerical issue
1074!||====================================================================
1075!|| zero2small ../starter/source/materials/tools/nlsqf.F
1076!||--- called by ------------------------------------------------------
1077!|| law69_guess1 ../starter/source/materials/tools/nlsqf.F
1078!||====================================================================
1079 SUBROUTINE zero2small(value, small)
1080#include "implicit_f.inc"
1081 my_real value, small
1082
1083 if ( abs(value) < small ) then
1084 if (value >= 0.0) then
1085 value = small
1086 else
1087 value = -small
1088 endif
1089 endif
1090 END SUBROUTINE zero2small
1091
1092!> @purpose
1093!> get one guess for starting Levenberg-Marquardt method
1094!> @param[in] LAWID material type(1=Ogden, 2=Mooney-Rivlin)
1095!> @param[out] NGUESS Number of returned guess
1096!> 0 = no return,already go through all combinations,
1097!> 1 = return one guess point successfully
1098!> @param[out] A0(NMUAL) returned guess
1099!||====================================================================
1100!|| law69_guess1 ../starter/source/materials/tools/nlsqf.F
1101!||--- called by ------------------------------------------------------
1102!|| law69_nlsqf ../starter/source/materials/tools/nlsqf.F
1103!|| law69_nlsqf_auto ../starter/source/materials/tools/law69_nlsqf_auto.F
1104!||--- calls -----------------------------------------------------
1105!|| law69_check ../starter/source/materials/tools/nlsqf.F
1106!|| rand1 ../starter/source/materials/tools/nlsqf.F
1107!|| zero2small ../starter/source/materials/tools/nlsqf.F
1108!||====================================================================
1109 SUBROUTINE law69_guess1(
1110 $ LAWID, NMUAL, ICHECK, MU_INCR, IRND1,
1111 $ A0, IA, MCOF_MIN, MCOF_MAX, NGUESS)
1112#include "implicit_f.inc"
1113C-----------------------------------------------
1114C A n a l y s e M o d u l e
1115C-----------------------------------------------
1116
1117 REAL RAND1
1118 EXTERNAL RAND1
1119
1120 INTEGER LAWID, NMUAL, ICHECK, MU_INCR, IRND1, NGUESS
1121 my_real A0(NMUAL)
1122 INTEGER IA(NMUAL)
1123
1124 my_real small
1125 data small /1e-15/
1126 save small
1127
1128 my_real mcof_min(nmual) ! lower bounds of material parameters
1129 my_real mcof_max(nmual) ! upper bounds of material parameters
1130
1131 real*4 randnum
1132
1133 INTEGER I, IVALID
1134
1135 nguess = 0
1136
1137 IF (lawid == 2) THEN
1138 a0(2) = 2.0
1139 a0(4) = -2.0
1140 ENDIF
1141
1142 DO WHILE (.true.)
1143 DO i = 1, nmual
1144 IF (ia(i) > 0) THEN
1145 randnum = rand1(irnd1)
1146 a0(i) = (1.0-randnum)*mcof_min(i) + randnum*mcof_max(i)
1147 CALL zero2small(a0(i), small)
1148 ENDIF
1149 ENDDO
1150
1151 CALL law69_check( lawid, a0, nmual,
1152 $ icheck, mu_incr, ivalid)
1153 IF (ivalid == 1) THEN
1154 nguess = 1
1155 GOTO 299
1156 ENDIF
1157 ENDDO
1158
1159299 CONTINUE
1160
1161 END SUBROUTINE law69_guess1
1162
1163!> @purpose
1164!> check if the guess of material properties is valid for LM optimization
1165!> @param[in] ICHECK checking level
1166!> <=0 no constraint (IVALID always return as 1)
1167!> 1 SUM( mu(i) * alpha(i) ) > 0.0
1168!> 2 mu(i) * alpha(i) > 0.0
1169!> @param[in] MU_INCR if condition mu(i) <= mu(i+1) is checked
1170!> 0=no, 1=yes
1171!> @param[out] IVALID 0=not valid, 1=valid
1172!||====================================================================
1173!|| law69_check ../starter/source/materials/tools/nlsqf.F
1174!||--- called by ------------------------------------------------------
1175!|| law69_guess1 ../starter/source/materials/tools/nlsqf.F
1176!|| law69_nlsqf ../starter/source/materials/tools/nlsqf.F
1177!|| law69_nlsqf_auto ../starter/source/materials/tools/law69_nlsqf_auto.F
1178!||====================================================================
1179 SUBROUTINE law69_check(LAWID, A, NMUAL,
1180 $ ICHECK, MU_INCR, IVALID)
1181#include "implicit_f.inc"
1182#include "units_c.inc"
1183 INTEGER LAWID, NMUAL, ICHECK, MU_INCR, IVALID
1184 my_real A(NMUAL)
1185
1186c local vars
1187 INTEGER J
1188 my_real SUM_MU_X_AF
1189
1190 IF (IDEBUG > 0) THEN
1191 WRITE(iout,*) 'LAW69_CHECK ...'
1192 DO j = 1, nmual
1193 WRITE(iout,*) 'J, A(J) = ', j, a(j)
1194 ENDDO
1195 ENDIF
1196
1197 ivalid = 1
1198 IF (icheck == 0) GOTO 49
1199
1200 IF (mu_incr > 0) THEN
1201 IF (lawid == 1) THEN
1202 IF (nmual >= 4) THEN
1203 DO j = 1, nmual/2-1
1204 IF ( a(2*j-1) > a(2*(j+1)-1) ) THEN
1205 IF (idebug > 0) THEN
1206 WRITE(iout,*) 'J,2*J-1,2*(J+1)-1 = ',
1207 $ j,2*j-1,2*(j+1)-1
1208 WRITE(iout,*) ' A(2*J-1), A(2*(J+1)-1) = ',
1209 $ a(2*j-1), a(2*(j+1)-1)
1210 WRITE(iout,*) ' MU(i) > MU(i+1) '
1211 ENDIF
1212 ivalid = 0
1213 goto 49
1214 ENDIF
1215 ENDDO
1216 ENDIF
1217 ENDIF
1218 ENDIF
1219
1220 IF (icheck == 1) THEN
1221C CHECK IF SUM(MU_I * ALPHA_I) > 0.0
1222 sum_mu_x_af = 0.0
1223 DO j = 1, nmual/2
1224 sum_mu_x_af = sum_mu_x_af + a(2*j-1) * a(2*j)
1225 ENDDO
1226 IF ( sum_mu_x_af <= 0.0) THEN
1227 IF (idebug > 0) THEN
1228 WRITE(iout,*) ' SUM_MU_X_AF = ', sum_mu_x_af
1229 ENDIF
1230 ivalid = 0
1231 GOTO 49
1232 ENDIF
1233 ELSEIF (icheck == 2) THEN
1234 DO j = 1, nmual/2
1235 IF ( ( a(2*j-1) * a(2*j) ) <= 0.0 ) THEN
1236 IF (idebug > 0) THEN
1237 WRITE(iout,*) ' A(2*J-1) * A(2*J) <= 0 '
1238 WRITE(iout,*) ' J, A(2*J-1), A(2*J) = ',
1239 $ j, a(2*j-1), a(2*j)
1240 ENDIF
1241 ivalid = 0
1242 GOTO 49
1243 ENDIF
1244 ENDDO
1245 ENDIF
1246
124749 CONTINUE
1248
1249 IF (idebug > 0) THEN
1250 WRITE(iout,*) 'IVALID = ', ivalid
1251 ENDIF
1252
1253 RETURN
1254
1255 END SUBROUTINE law69_check
1256
1257!> @purpose
1258!> Esstimate lwer/upper bounds of mu_i and alpha_i
1259!> which will be used to generate initial guess
1260!> @param[in] LAWID 1=Ogden, 2=Mooney-Rivlin
1261!> @param[in] ICHECK Checking level
1262!> @param[in] X(NPT) stretch values of given stress/strain(stretch) curve
1263!> @param[in] Y(NPT) engineering stress values of given stress/strain curve
1264!> @param[in] NPT number of points on the given stress/strain curve
1265
1266!> @param[out] MCOF_MIN(NMUAL) estimated lower bounds of A(1:NMUAL)
1267!> @param[out] MCOF_MAX(NMUAL) estimated upper bounds of A(1:NMUAL)
1268!||====================================================================
1269!|| law69_guess_bounds ../starter/source/materials/tools/nlsqf.F
1270!||--- called by ------------------------------------------------------
1271!|| law69_nlsqf ../starter/source/materials/tools/nlsqf.F
1272!|| law69_nlsqf_auto ../starter/source/materials/tools/law69_nlsqf_auto.F
1273!||====================================================================
1274 SUBROUTINE law69_guess_bounds(LAWID, ICHECK, X, Y, NPT,
1275 $ NMUAL, MCOF_MIN, MCOF_MAX,ICOMP)
1276#include "implicit_f.inc"
1277#include "units_c.inc"
1278
1279 INTEGER LAWID, ICHECK, NPT,ICOMP
1280 my_real X(NPT), Y(NPT)
1281 INTEGER NMUAL
1282 my_real MCOF_MIN(NMUAL), MCOF_MAX(NMUAL)
1283
1284 INTEGER I
1285
1286 my_real AVE_SLOPE, mu_max, dx
1287
1288 ave_slope = zero
1289 DO i = 1, npt
1290 dx = x(i) - one
1291c avolid dx to be too small
1292 IF (dx >= zero) THEN
1293 dx = max(dx, em6)
1294 ELSE
1295 dx = min(dx, -em6)
1296 ENDIF
1297 ave_slope = ave_slope + (y(i) - zero) / dx
1298!ok AVE_SLOPE = MAX(AVE_SLOPE , (Y(I) - ZERO) / dx )
1299 ENDDO
1300 ave_slope = ave_slope / (one * npt)
1301 IF (idebug > 0) THEN
1302 WRITE(iout, *) 'AVE_SLOPE = ', ave_slope
1303 ENDIF
1304 mu_max = ave_slope
1305 IF(icomp == 0 ) mu_max = max(ave_slope, twenty)
1306 IF (lawid == 1) THEN
1307 DO i = 1, nmual/2
1308 ! mu_i
1309 mcof_min(2*i-1) = -mu_max
1310 mcof_max(2*i-1) = mu_max
1311
1312 ! alpha_i
1313 mcof_min(2*i) = -10.0
1314 mcof_max(2*i) = +10.0
1315 ENDDO
1316 ELSEIF (lawid == 2) THEN
1317
1318 IF (icheck == 2) THEN
1319 mcof_min(1) = 0.0
1320 mcof_max(1) = mu_max
1321 ELSE
1322 mcof_min(1) = -mu_max
1323 mcof_max(1) = mu_max
1324 ENDIF
1325
1326 mcof_min(2) = +2.0
1327 mcof_max(2) = +2.0
1328
1329 IF (icheck == 2) THEN
1330 mcof_min(3) = -mu_max
1331 mcof_max(3) = 0.0
1332 ELSE
1333 mcof_min(3) = -mu_max
1334 mcof_max(3) = mu_max
1335 ENDIF
1336
1337 mcof_min(4) = -2.0
1338 mcof_max(4) = -2.0
1339
1340 ENDIF
1341
1342 IF (idebug > 0) THEN
1343 DO i = 1, nmual
1344 WRITE(iout,*) 'I, MCOF_MIN(I), MCOF_MAX(I) = ',
1345 $ i, mcof_min(i), mcof_max(i)
1346 ENDDO
1347 ENDIF
1348
1349 END SUBROUTINE law69_guess_bounds
1350
1351!> @purpose
1352!> return a random value within [0.0, 1.0) with uniform distribution
1353!> @param[in] IVAL integer for calculating next random number
1354!> @note
1355!> - random number generator of Park and Miller.
1356!> - Schrage's modification is used to avoid 'a*(m-1)' overflows for 32-bit.
1357!||====================================================================
1358!|| rand1 ../starter/source/materials/tools/nlsqf.F
1359!||--- called by ------------------------------------------------------
1360!|| law69_guess1 ../starter/source/materials/tools/nlsqf.F
1361!||--- calls -----------------------------------------------------
1362!|| my_exit ../starter/source/output/analyse/analyse.c
1363!||====================================================================
1364 REAL function rand1(ival)
1365#include "implicit_f.inc"
1366 INTEGER ival
1367 INTEGER ia, im, iq, ir
1368 REAL rim
1369 PARAMETER (ia=69621, im=2147483647,rim=1./im,iq=30845,ir=23902)
1370
1371 INTEGER idq
1372
1373 IF (ival <= 0) THEN
1374 WRITE(*,*) 'ERROR, IVAL <= 0'
1375 WRITE(*,*) __file__,__line__
1376 CALL my_exit(2)
1377 ELSEIF (ival > im) THEN
1378 WRITE(*,*) 'ERROR, IVAL > IM'
1379 WRITE(*,*) __file__,__line__
1380 CALL my_exit(2)
1381 ENDIF
1382
1383 idq = ival / iq
1384 ival = ia * (ival - idq * iq) - ir * idq
1385 if (ival < 0) ival = ival + im
1386
1387 rand1 = real(ival) * rim
1388
1389 RETURN
1390 END FUNCTION rand1
my_exit
void my_exit(int *i)
Definition analyse.c:1038
my_real
#define my_real
Definition cppsort.cpp:32
alpha
#define alpha
Definition eval.h:35
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
message_mod
Definition message_mod.F:1257
names_and_titles_mod
Definition names_and_titles_mod.F:1002
names_and_titles_mod::nchartitle
integer, parameter nchartitle
Definition names_and_titles_mod.F:1008
zero2small
subroutine zero2small(value, small)
Definition nlsqf.F:1080
ogden1
subroutine ogden1(x, a, dyda, ia, na)
@purpose calculate d(stress)/d(A) OGDEN (uniaxial tension)
Definition nlsqf.F:508
rand1
real function rand1(ival)
@purpose return a random value within [0.0, 1.0) with uniform distribution
Definition nlsqf.F:1365
mrqmin
subroutine mrqmin(x, y, sig, ndata, a, ia, ma, covar, alpha, nca, errnow, ifuncs, gamma, iret, icomp, mmax, atry)
Definition nlsqf.F:45
law69_nlsqf
subroutine law69_nlsqf(lawid, stretch, y, ma, nmual, npt, mual, icheck, nstart, errtol, id, titr, is_encrypted)
Definition nlsqf.F:565
inversion
subroutine inversion(a, n, np, b, m, mp, iret)
Definition nlsqf.F:304
law69_check
subroutine law69_check(lawid, a, nmual, icheck, mu_incr, ivalid)
@purpose check if the guess of material properties is valid for LM optimization
Definition nlsqf.F:1181
law69_guess_bounds
subroutine law69_guess_bounds(lawid, icheck, x, y, npt, nmual, mcof_min, mcof_max, icomp)
@purpose Esstimate lwer/upper bounds of mu_i and alpha_i which will be used to generate initial guess
Definition nlsqf.F:1276
mrqcof
subroutine mrqcof(x, y, sig, ndata, a, ia, ma, alpha, beta, nalp, errnow, ifuncs, icomp)
Definition nlsqf.F:188
covsrt
subroutine covsrt(covar, sizcovar, ma, ia, mfit)
Definition nlsqf.F:420
law69_guess1
subroutine law69_guess1(lawid, nmual, icheck, mu_incr, irnd1, a0, ia, mcof_min, mcof_max, nguess)
@purpose get one guess for starting Levenberg-Marquardt method
Definition nlsqf.F:1112
ogden0
subroutine ogden0(x, a, y, na)
@purpose calculate stress from stretch for OGDEN (uniaxial tension)
Definition nlsqf.F:473
ancmsg
subroutine ancmsg(msgid, msgtype, anmode, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16, i17, i18, i19, i20, r1, r2, r3, r4, r5, r6, r7, r8, r9, c1, c2, c3, c4, c5, c6, c7, c8, c9, prmode)
Definition message.F:895