1!
2!     CalculiX - A 3-dimensional finite element program
3!              Copyright (C) 1998-2021 Guido Dhondt
4!
5!     This program is free software; you can redistribute it and/or
6!     modify it under the terms of the GNU General Public License as
7!     published by the Free Software Foundation(version 2);
8!
9!
10!     This program is distributed in the hope that it will be useful,
11!     but WITHOUT ANY WARRANTY; without even the implied warranty of
12!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13!     GNU General Public License for more details.
14!
15!     You should have received a copy of the GNU General Public License
16!     along with this program; if not, write to the Free Software
17!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
18!
19      subroutine objective_modalstress(ndesi,neq,b,daldx,bfix,jqs,
20     &     irows,df,iev,nev,z,dgduz,d,iobject,dgdx,dfm)
21!
22      implicit none
23!
24      integer i,j,idesvar,ndesi,neq(*),irows(*),iev,nev,iobject,jqs(*)
25!
26      real*8 b(*),daldx(*),bfix(*),df(*),sensi,sum,z(neq(2),*),
27     &     d(*),dgdx(ndesi,*),dgduz(*),dfm(*)
28!
29!     calculation of the modal stress sensitivity
30!
31      do idesvar=1,ndesi
32!
33!       dlambda_i/ds*M*U_i
34!
35        do j=1,neq(2)
36          b(j)=daldx(idesvar)*bfix(j)
37        enddo
38!
39!       dF/ds := (-dK/ds+lambda_i*dM/ds+dlambda_i/ds*M)*U_i
40!
41        do j=jqs(idesvar),jqs(idesvar+1)-1
42          b(irows(j))=b(irows(j))+df(j)
43        enddo
44!
45        sensi=0.d0
46        do i=1,nev
47          if(i.eq.iev+1) then
48!
49!           calculation of U_i^T*(dM/ds*U_i) = -2*betas
50!
51            sum=0.d0
52            do j=jqs(idesvar),jqs(idesvar+1)-1
53              sum=sum+z(irows(j),i)*dfm(j)
54            enddo
55            sensi=sensi-sum*dgduz(i)/2.d0
56            cycle
57          endif
58!
59!         calculation of U_i^T . dF/ds
60!
61          sum=0.d0
62          do j=1,neq(2)
63            sum=sum+z(j,i)*b(j)
64          enddo
65          sensi=sensi+sum*dgduz(i)/(d(iev+1)-d(i))
66        enddo
67        dgdx(idesvar,iobject)=dgdx(idesvar,iobject)+sensi
68      enddo
69!
70      return
71      end
72