Import of IFEM source files from the old CVS repository on afem

git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@769 e10b68d5-8a6e-419e-a041-bce267b0401d
This commit is contained in:
akva
2011-02-08 13:48:45 +00:00
committed by Knut Morten Okstad
parent ad9c648c8f
commit c88a28533d
214 changed files with 58868 additions and 0 deletions

View File

@@ -0,0 +1,71 @@
# $Id: Makefile,v 1.7 2011-02-08 09:32:18 kmo Exp $
#===============================================================================
#
# File: Makefile
#
# Created: Jan 12 2010
#
# Author: Knut Morten Okstad / SINTEF
#
# Description:
# Makefile for spline-based isogeometric finite element simulation applications.
#
#===============================================================================
CSRC = main_NonLinEl.C NonlinearElasticityTL.C NonlinearElasticityUL.C \
NonlinearElasticityULMX.C NonlinearElasticityULMixed.C \
NonlinearElasticity.C NeoHookeElasticity.C \
SIMFiniteDefEl.C ../SIMLinEl2D.C ../SIMLinEl3D.C
FSRC = accKM.f accKMx.f stiff_TL.f
OBJS = $(CSRC:.C=.o) $(FSRC:.f=.o)
APP = nonlinel
CXX = g++
FC = gfortran
MYLIB = ../../libSplineFEM.a libMaterial.a
GOTLIB = -L../../GoTools/lib -lGoTrivariate -lGoToolsCore
VTFLIB = -L../../GLviewExpressWriter/lib -lVTFExpressAPI
BLASLIB = -larpack -llapack -lblas
SLULIB = -lsuperlu
#SLULIB = -L../../SuperLU_MT/lib -lsuperlu_mt_PTHREAD -lpthread
#SAMGLIB = -L../../SAMG/samg -lamg_coo
LIBS = $(MYLIB) $(GOTLIB) $(VTFLIB) $(SAMGLIB) $(SLULIB) $(BLASLIB)
#CDEBUG = -g -DINDEX_CHECK=2 -DINT_DEBUG=2
CDEBUG = -O2
FFLAGS = -O3 -funroll-loops
CXXFLAGS = $(CDEBUG) -Dreal=double -DUSE_CBLAS -DUSE_FTNMAT \
-I. -I.. -I../../Integrands -I../../SIM -I../../ASM \
-I../../LinAlg -I../../Utility
all: $(APP)
$(APP): $(OBJS) $(MYLIB) Makefile
$(CXX) -o $@ $(OBJS) $(LIBS) -lgfortran
../../libSplineFEM.a:
cd ../..; make
libMaterial.a:
cd Material; make
clean:
rm -f $(OBJS)
sed '1,/^# DO NOT DELETE/w Makefile.tmp' < Makefile > /dev/null
mv Makefile.tmp Makefile
depend:
sed '1,/^# DO NOT DELETE/w Makefile.tmp' < Makefile > /dev/null
$(CXX) -M $(CXXFLAGS) *.C >> Makefile.tmp
mv Makefile Makefile.bak
mv Makefile.tmp Makefile
tar: clean
tar cfv FDSolve.tar Makefile $(FSRC) $(CSRC) *.h \
Material/Makefile Material/*.[fh]
gzip FDSolve.tar
# DO NOT DELETE

View File

@@ -0,0 +1,23 @@
# $Id: Makefile,v 1.1 2010-12-18 20:03:54 kmo Exp $
LIB = ../libMaterial.a
SRC = amat3d.f cons2d.f cons3d.f eigs3d.f ieps3d.f liel3d.f mdma3d.f \
monh3d.f pcst3d.f push3d.f stnh3d.f vfnh3d.f wlog3d.f rprin0.f
OBJS = $(SRC:.f=.o)
FC = gfortran
FFLAGS = -O3
#FFLAGS = -g
all: feninc $(LIB)
clean:
-rm $(OBJS)
$(LIB): $(OBJS)
$(AR) $(ARFLAGS) $@ $(OBJS)
feninc:
ln -s . include
ln -s . feninc

View File

@@ -0,0 +1,162 @@
subroutine amat3d (taup, aap, q, tau, aa,
& ipswb3, iwr, ierr)
C
C ---------------------------------------------------------------------
C
C F E N R I S ROUTINE: AMAT3D NTNU/SINTEF/DNV
C
C PURPOSE:
C AMAT3D: Transform Kirchhoff stresses and material moduli from
C principal to standard basis
C
C METHOD:
C
C
C ARGUMENTS INPUT:
C taup(3) - Principal deviatoric Kirchhoff stresses
C aap(6,6) - Deviatoric tangent moduli in principal basis
C q(6,6) - Transformation matrix: [principal]-->[std] basis
C
C ARGUMENTS OUTPUT:
C tau(3) - Deviatoric Kirchhoff stresses (ref to standard basis)
C aa(6,6) - Tangent moduli in standard basis
C ierr - Error flag
C = 1 :
C
C COMMON BLOCKS:
C const
C
C PERIPHERAL UNITS:
C None
C
C INTERNAL VARIABLES:
C
C
C
C PRINT SWITCH:
C ipswb3 = 0 Gives no print
C ipswb3 = 2 Gives enter and leave
C ipswb3 = 3 Gives in addition parameters on input
C ipswb3 = 5 Gives in addition parameters on output
C
C LIMITS:
C
C CODING:
C STANDART ADHERED TO: Computas report no 78-949
C AUTHOR: Kjell Magne Mathisen
C NTNU, Department of structural engineering
C November 14, 2010
C
C INSPECTED:
C
C REVISIONS:
C xx-OCT-2010 K.M.Mathisen NTNU
C Changed XG1-XG8 to XG and UG1-UG8 to UG.
C
C ---------------------------------------------------------------------
C
implicit none
C
integer i, j, k, ipswb3, iwr, ierr
C
real*8 taup(3), aap(6,6), q(6,6), tau(6), aa(6,6), vi(6)
C
include 'include/feninc/const.h'
C
C Entry section
C
ierr = 0
C
if (ipswb3 .gt. 0) then
write(iwr,9010) 'ENTERING SUBROUTINE AMAT3D'
if (ipswb3 .gt. 2) then
write(iwr,9010) 'WITH INPUT ARGUMENTS'
call rprin0(taup, 1, 3, 'taup ', iwr)
call rprin0(aap , 6, 6, 'aap ', iwr)
call rprin0(q , 6, 6, 'q ', iwr)
endif
endif
C
C Compute stress tensor
C
do j = 1,6
tau(j) = q(j,1)*taup(1) + q(j,2)*taup(2) + q(j,3)*taup(3)
end do ! j
C
C Compute tangent matrix in standard basis
C
do i = 1,6
C
C Left multiplication: V = Q * AL
C
do k = 1,3
vi(k) = zero
do j = 1,3
vi(k) = vi(k) + q(i,j)*aap(j,k)
end do ! j
end do ! k
C
do k = 4,6
vi(k) = q(i,k)*aap(k,k)
end do ! k
C
C Right multiplication: AA = V * Q_t
C
do j = 1,i
aa(i,j) = zero
do k = 1,6
aa(i,j) = aa(i,j) + vi(k)*q(j,k)
end do ! k
C
C Fill in symmetric part
C
aa(j,i) = aa(i,j)
C
end do ! j
C
end do ! i
go to 8000
C
C Error section
C
7000 continue
ierr =-1
go to 7900
7100 continue
ierr =-2
go to 7900
C ----------
7900 continue
write(iwr,9020) '*** ERROR IN SUBROUTINE AMAT3D *** IERR = ', ierr
write(iwr,9010) 'WITH INPUT ARGUMENTS'
call rprin0(taup, 1, 3, 'taup ', iwr)
call rprin0(aap , 6, 6, 'aap ', iwr)
call rprin0(q , 6, 6, 'q ', iwr)
C
write(iwr,9010) 'WITH OUTPUT ARGUMENTS'
call rprin0(tau, 1, 3, 'tau ', iwr)
call rprin0(aa , 6, 6, 'aa ', iwr)
go to 8010
C
C Closing section
C
8000 continue
C
if (ipswb3 .gt. 0) then
write(iwr,9010) 'LEAVING SUBROUTINE AMAT3D'
if (ipswb3 .gt. 3) then
write(iwr,9010) 'WITH OUTPUT ARGUMENTS'
call rprin0(tau, 1, 3, 'tau ', iwr)
call rprin0(aa , 6, 6, 'aa ', iwr)
endif
endif
C
8010 continue
C
return
C
9010 format(3X,A)
9020 format(3X,A,I6)
9030 format(3X,A,E15.6)
C
end

View File

@@ -0,0 +1,29 @@
subroutine cons2d (mTYP,mVER,detF,F,pmat,U,Sig,Cst,ipsw,iwr,ierr)
C
C CONS2D: 2D wrapper on the 3D material of FENRIS.
C
implicit none
C
integer mTYP, mVER, ipsw, iwr, ierr
real*8 detF, U, F(4), pmat(*), Sig(3), Cst(3,3)
real*8 F3D(3,3), S3D(6), C3D(6,6)
C
F3D = 0.0D0
F3D(1,1) = F(1)
F3D(2,1) = F(2)
F3D(1,2) = F(3)
F3D(2,2) = F(4)
F3D(3,3) = 1.0D0
call cons3d (mTYP,mVER,detF,F3D,pmat,U,S3D,C3D,ipsw,iwr,ierr)
Cst(1:2,1:2) = C3D(1:2,1:2)
Cst(1:2,3) = C3D(1:2,4)
Cst(3,1:2) = C3D(4,1:2)
Cst(3,3) = C3D(4,4)
Sig(1:2) = S3D(1:2)
Sig(3) = S3D(4)
C
return
end

View File

@@ -0,0 +1,226 @@
subroutine cons3d (mTYP, mVER, detF, F, pmat, Engy, Sig, Cst,
& ipswb3, iwr, ierr)
C
C ---------------------------------------------------------------------
C
C F E N R I S ROUTINE: CONS3D NTNU/SINTEF/DNV
C
C PURPOSE:
C CONS3D: Driver for finite deformation constitutive models for
C 3D solids
C
C METHOD:
C Compute Cauchy stresses and spatial constitutive tensor
C for finite deformation constitutive models for 3D solids.
C
C The following material models are available:
C mTYP = 10 Linear elastic material, 3-D stress
C mTYP = 20 Hyperelastic compressible neoHookean model, 3-D stress
C mTYP = 40 Elastic-plastic strain hardening material,
C combined isotropic and kinematic hardening,
C 3-D stress
C
C ARGUMENTS INPUT:
C mTYP - Material type number
C mVER - Material version number
C detF - Determinant of the deformation gradient
C F - Deformation gradient
C pmat(2) - Material property array
C
C ARGUMENTS OUTPUT:
C Engy - Strain energy density
C Sig(6) - Cauchy stresses (spatial stresses)
C Cst(6,6) - Spatial constitutive tensor
C ierr - Error flag
C = -1 : Illegal number of element nodes
C
C COMMON BLOCKS:
C const
C
C PERIPHERAL UNITS:
C None
C
C INTERNAL VARIABLES:
C Bmod - Material property (Bulk modulus)
C Emod - Material property (Youngs modulus)
C Pratio - Material property (Poissons ratio)
C Smod - Material property (Shear modulus)
C incomp - logical variable
C = .true. if incompressible material ): Pratio = 0
C
C PRINT SWITCH:
C ipswb3 = 0 Gives no print
C ipswb3 = 2 Gives enter and leave
C ipswb3 = 3 Gives in addition parameters on input
C ipswb3 = 5 Gives in addition parameters on output
C
C LIMITS:
C
C CODING:
C STANDART ADHERED TO: Computas report no 78-949
C AUTHOR: Kjell Magne Mathisen
C NTNU, Department of structural engineering
C September 14, 2010
C
C INSPECTED:
C
C REVISIONS:
C xx-OCT-2010 K.M.Mathisen NTNU
C Changed XG1-XG8 to XG and UG1-UG8 to UG.
C
C ---------------------------------------------------------------------
C
implicit none
C
integer mTYP, mVER, ipswb3, iwr, ierr
real*8 detF, Engy, F(9), pmat(2), Sig(6), Cst(6,6)
real*8 Bmod, Emod, Pratio, Smod
logical incomp
C
include 'include/feninc/const.h'
C
C Entry section
C
ierr = 0
C
if (ipswb3 .gt. 0) then
write(iwr,9010) 'ENTERING SUBROUTINE CONS3D'
if (ipswb3 .gt. 2) then
write(iwr,9010) 'WITH INPUT ARGUMENTS'
write(iwr,9020) 'mTYP =', mTYP
write(iwr,9020) 'mVER =', mVER
write(iwr,9030) 'detF =', detF
write(iwr,9030) 'Pmat =', pmat
call rprin0(F , 3, 3, 'F ', iwr)
endif
endif
C
C
C Compute Shear and Bulk modulus
C
if (mTYP .lt. 0) then
Bmod = pmat(1)
Smod = pmat(2)
else
Emod = pmat(1)
Pratio = pmat(2)
Smod = one2 * Emod / (one + Pratio)
incomp = Pratio .eq. one2
if (incomp) then
C
C Incompressible material
C
Bmod = zero
else
Bmod = one3 * Emod / (one - two * Pratio)
end if
end if
C
C Compute Cauchy stresses and spatial constitutive tensor
C
if (abs(mTYP) .eq. 10) then
C
if (mVER .eq. 0) then
C
C Linear elastic isotropic material
C
call liel3d(detF, F, Bmod, Smod, Engy, Sig, Cst,
& ipswb3, iwr, ierr)
if (ierr .lt. 0) go to 7000
C
elseif (mVER .gt. 0 .and. mVER .lt. 5) then
C
C Standard hyperelastic neoHookean model
C
C Compute Lame's parameter for compressibel material
C
if (.not. incomp) Bmod = Bmod - two3 * Smod
C
call stnh3d(mVER, detF, F, Bmod, Smod, Engy, Sig, Cst,
& ipswb3, iwr, ierr)
if (ierr .lt. 0) go to 7000
C
elseif (mVER .gt. 10 .and. mVER .lt. 15) then
C
C Modified hyperelastic neoHookean model
C
call monh3d(mVER-10, detF, F, Bmod, Smod, Engy, Sig, Cst,
& ipswb3, iwr, ierr)
if (ierr .lt. 0) go to 7000
C
elseif (mVER .eq. 20) then
C
C Isotropic elasticity formulated in terms of
C principal logarithmic stretches
C
call ieps3d(detF, F, Bmod, Smod, Engy, Sig, Cst,
& ipswb3, iwr, ierr)
if (ierr .lt. 0) go to 7000
C
else
C
C Illegal material version
C
go to 7100
endif
C
else
C
C Illegal material type
C
go to 7200
endif
go to 8000
C
C Error section
C
7000 continue
ierr =-1
go to 7900
7100 continue
ierr =-2
go to 7900
7200 continue
ierr =-3
go to 7900
C ----------
7900 continue
write(iwr,9020) '*** ERROR IN SUBROUTINE CONS3D *** IERR = ', ierr
write(iwr,9010) 'WITH INPUT ARGUMENTS'
write(iwr,9020) 'mTYP =', mTYP
write(iwr,9020) 'mVER =', mVER
write(iwr,9030) 'detF =', detF
write(iwr,9030) 'Pmat =', pmat
call rprin0(F , 3, 3, 'F ', iwr)
C
write(iwr,9010) 'WITH OUTPUT ARGUMENTS'
write(iwr,9030) 'Engy =', Engy
call rprin0(Sig , 1, 6, 'Sig ', iwr)
call rprin0(Cst , 6, 6, 'Cst ', iwr)
call flush(iwr)
go to 8010
C
C Closing section
C
8000 continue
C
if (ipswb3 .gt. 0) then
write(iwr,9010) 'LEAVING SUBROUTINE CONS3D'
if (ipswb3 .gt. 3) then
write(iwr,9010) 'WITH OUTPUT ARGUMENTS'
write(iwr,9030) 'Engy =', Engy
call rprin0(Sig , 1, 6, 'Sig ', iwr)
call rprin0(Cst , 6, 6, 'Cst ', iwr)
endif
call flush(iwr)
endif
C
8010 continue
C
return
C
9010 format(3X,A)
9020 format(3X,A,I6)
9030 format(3X,A,1P,2E15.6)
C
end

View File

@@ -0,0 +1,6 @@
double precision zero,one,two,four
double precision one2,one3,one4,one6,onep5,two3,four3
parameter ( zero=0.0D0, one=1.0D0, two=2.0D0, four=4.0D0 )
parameter ( one2=0.5D0, one3=one/3.0D0, one4=0.25D0 )
parameter ( one6=one/6.0D0, onep5=1.5D0 )
parameter ( two3=two/3.0D0, four3=four/3.0D0 )

View File

@@ -0,0 +1,257 @@
subroutine eigs3d (v, d, nrotd,
& ipswb3, iwr, ierr)
C
C ---------------------------------------------------------------------
C
C F E N R I S ROUTINE: EIGS3D NTNU/SINTEF/DNV
C
C PURPOSE:
C EIGS3D: Compute eigenvalues/vectors for 3x3 symmetric matrix
C
C METHOD:
C Storage done as follows:
C
C | v(1,1) v(1,2) v(1,3) | | d(1) a(1) a(3) |
C | v(2,1) v(2,2) v(2,3) | = | a(1) d(2) a(2) |
C | v(3,1) v(3,2) v(3,3) | | a(3) a(2) d(3) |
C
C Transformations performed on d(i) and a(i) and v(i,j) become
C eigenvectors. Thus, original array is destroyed.
C
C ARGUMENTS INPUT:
C v(3,3) - Matrix for which eigenvalues/vectors should be
C computed for
C
C ARGUMENTS OUTPUT:
C v(3,3) - Matrix of eigenvectors (by column)
C d(3) - Eigenvalues associated with the columns of v
C nrotd - number of rotations to diagonalize
C ierr - Error flag
C = -1 :
C
C COMMON BLOCKS:
C const
C
C PERIPHERAL UNITS:
C None
C
C INTERNAL VARIABLES:
C
C
C
C PRINT SWITCH:
C ipswb3 = 0 Gives no print
C ipswb3 = 2 Gives enter and leave
C ipswb3 = 3 Gives in addition parameters on input
C ipswb3 = 5 Gives in addition parameters on output
C
C LIMITS:
C
C CODING:
C STANDART ADHERED TO: Computas report no 78-949
C AUTHOR: Kjell Magne Mathisen
C NTNU, Department of structural engineering
C November 14, 2010
C
C INSPECTED:
C
C REVISIONS:
C xx-OCT-2010 K.M.Mathisen NTNU
C Changed XG1-XG8 to XG and UG1-UG8 to UG.
C
C ---------------------------------------------------------------------
C
implicit none
C
integer nrotd, i, its, j, k, ipswb3, iwr, ierr
C
real*8 g, h, aij, sm, thresh, t, c, s, tau,
* v(3,3), d(3), a(3), b(3), z(3)
C
include 'include/feninc/const.h'
C
C Entry section
C
ierr = 0
C
if (ipswb3 .gt. 0) then
write(iwr,9010) 'ENTERING SUBROUTINE EIGS3D'
if (ipswb3 .gt. 2) then
write(iwr,9010) 'WITH INPUT ARGUMENTS'
call rprin0(v, 3, 3, 'v ', iwr)
endif
endif
C
C Move array into one-d arrays
C
a(1) = v(1,2)
a(2) = v(2,3)
a(3) = v(1,3)
C
do i = 1,3
C
d(i) = v(i,i)
b(i) = d(i)
z(i) = zero
C
do j = 1,3
C
v(i,j) = zero
C
end do ! j
C
v(i,i) = one
C
end do ! i
C
C Check for diagonal case
C
sm = abs(a(1)) + abs(a(2)) + abs(a(3))
g = abs(d(1)) + abs(d(2)) + abs(d(3))
C
if (sm .lt. 1.d-13*g) go to 8000
C
nrotd = 0
C
do its = 1, 50
C
C Set convergence test and threshold
C
sm = abs(a(1)) + abs(a(2)) + abs(a(3))
if (sm .eq. zero) go to 8000
C
if (its .lt. 4) then
C
thresh = 0.011d0*sm
C
else
C
thresh = zero
C
end if
C
C Perform sweeps for rotations
C
do i = 1,3
C
j = mod(i,3) + 1
k = mod(j,3) + 1
C
aij = a(i)
g = 100.d0*abs(aij)
C
if (abs(d(i))+g .ne. abs(d(i)) .or.
& abs(d(j))+g .ne. abs(d(j))) then
C
if (abs(aij) .gt. thresh) then
C
a(i) = zero
h = d(j) - d(i)
C
if (abs(h)+g .eq. abs(h)) then
t = aij/h
else
t = sign(two,h/aij)/(abs(h/aij)
& + sqrt(four+(h/aij)**2))
end if
C
C Set rotation parameters
C
c = one/sqrt(one+t*t)
s = t*c
tau = s/(one+c)
C
C Rotate diagonal terms
C
h = t*aij
z(i) = z(i) - h
z(j) = z(j) + h
d(i) = d(i) - h
d(j) = d(j) + h
C
C Rotate off-diagonal terms
C
h = a(j)
g = a(k)
a(j) = h + s*(g - h*tau)
a(k) = g - s*(h + g*tau)
C
C Rotate eigenvectors
C
do k = 1,3
C
g = v(k,i)
h = v(k,j)
v(k,i) = g - s*(h + g*tau)
v(k,j) = h + s*(g - h*tau)
C
end do ! k
C
nrotd = nrotd + 1
C
end if
C
else
C
a(i) = zero
C
end if
C
end do ! i
C
C Update diagonal terms
C
do i = 1,3
C
b(i) = b(i) + z(i)
d(i) = b(i)
z(i) = zero
C
end do ! i
C
end do ! its
go to 8000
C
C Error section
C
7000 continue
ierr =-1
go to 7900
7100 continue
ierr =-2
go to 7900
C ----------
7900 continue
write(iwr,9020) '*** ERROR IN SUBROUTINE EIGS3D *** IERR = ', ierr
write(iwr,9010) 'WITH INPUT ARGUMENTS'
call rprin0(v, 3, 3, 'v ', iwr)
C
write(iwr,9010) 'WITH OUTPUT ARGUMENTS'
write(iwr,9020) 'nrotd =', nrotd
call rprin0(v, 3, 3, 'v ', iwr)
call rprin0(d, 1, 3, 'd ', iwr)
go to 8010
C
C Closing section
C
8000 continue
C
if (ipswb3 .gt. 0) then
write(iwr,9010) 'LEAVING SUBROUTINE EIGS3D'
if (ipswb3 .gt. 3) then
write(iwr,9010) 'WITH OUTPUT ARGUMENTS'
write(iwr,9020) 'nrotd =', nrotd
call rprin0(v, 3, 3, 'v ', iwr)
call rprin0(d, 1, 3, 'd ', iwr)
endif
endif
C
8010 continue
C
return
C
9010 format(3X,A)
9020 format(3X,A,I6)
9030 format(3X,A,E15.6)
C
end

View File

@@ -0,0 +1,247 @@
subroutine ieps3d (detF, F, lambda, mu, Engy, Sig, Cst,
& ipswb3, iwr, ierr)
C
C ---------------------------------------------------------------------
C
C F E N R I S ROUTINE: IEPS3D NTNU/SINTEF/DNV
C
C PURPOSE:
C IEPS3D: Compute Cauchy stresses and spatial constitutive tensor
C for a isotropic elastic 3D material formulated in
C principal stretches
C
C METHOD:
C
C
C
C ARGUMENTS INPUT:
C detF - Determinant of the deformation gradient
C F(3,3) - Deformation gradient
C lambda - Lame's constant (=K-2*G/3)
C mu - Shear modulus
C
C ARGUMENTS OUTPUT:
C Engy - Strain energy density
C Sig(6) - Cauchy stresses (spatial stresses)
C Cst(6,6) - Spatial constitutive tensor
C ierr - Error flag
C = 1 :
C
C COMMON BLOCKS:
C const
C
C PERIPHERAL UNITS:
C None
C
C INTERNAL VARIABLES:
C b(6) - Left Cauchy-Green deformation tensor
C
C
C PRINT SWITCH:
C ipswb3 = 0 Gives no print
C ipswb3 = 2 Gives enter and leave
C ipswb3 = 3 Gives in addition parameters on input
C ipswb3 = 5 Gives in addition parameters on output
C
C LIMITS:
C
C CODING:
C STANDART ADHERED TO: Computas report no 78-949
C AUTHOR: Kjell Magne Mathisen
C NTNU, Department of structural engineering
C November 14, 2010
C
C INSPECTED:
C
C REVISIONS:
C xx-OCT-2010 K.M.Mathisen NTNU
C Changed XG1-XG8 to XG and UG1-UG8 to UG.
C
C ---------------------------------------------------------------------
C
implicit none
C
integer ipswb3, iwr, ierr
integer i, j, k, l, iVF, i3, ierrl, nrotd
C
real*8 detF, Engy, F(3,3), lambda, mu, Sig(6), Cst(6,6)
real*8 aap(6,6), b(6), bpr(6), Cstv, detFi, EngU, EngW,
& epsd(3), Press, q(6,6), qen(3,3), taup(3)
C
include 'include/feninc/const.h'
C
C Entry section
C
ierr = 0
C
if (ipswb3 .gt. 0) then
write(iwr,9010) 'ENTERING SUBROUTINE IEPS3D'
if (ipswb3 .gt. 2) then
write(iwr,9010) 'WITH INPUT ARGUMENTS'
write(iwr,9030) 'detF =', detF
write(iwr,9030) 'lambda =', lambda
write(iwr,9030) 'mu =', mu
call rprin0(F , 3, 3, 'F ', iwr)
endif
endif
C
C Initialize Cauchy stresses and spatial constitutive tensor
C
Sig = zero
Cst = zero
C
C Compute Left Cauchy-Green deformation tensor
C
C b = F * F^T
C
b(1) = F(1,1)*F(1,1) + F(1,2)*F(1,2) + F(1,3)*F(1,3)
b(2) = F(2,1)*F(2,1) + F(2,2)*F(2,2) + F(2,3)*F(2,3)
b(3) = F(3,1)*F(3,1) + F(3,2)*F(3,2) + F(3,3)*F(3,3)
b(4) = F(1,1)*F(2,1) + F(1,2)*F(2,2) + F(1,3)*F(2,3)
b(5) = F(2,1)*F(3,1) + F(2,2)*F(3,2) + F(2,3)*F(3,3)
b(6) = F(1,1)*F(3,1) + F(1,2)*F(3,2) + F(1,3)*F(3,3)
C
C Compute principal stretches and directions
C
C Initialization: Upon entry to EIGS3D array qen(3,3) contains
C left Cauchy-Green tensor b; upon exit,
C qen(3,3) contains eigenvectors
C
do i = 1,3
C
i3 = i + 3
j = 1 + mod(i,3)
C
qen(i,i) = b(i)
qen(i,j) = b(i3)
qen(j,i) = b(i3)
C
end do ! i
C
call eigs3d(qen, bpr, nrotd, ipswb3, iwr, ierr)
if (ierrl .lt. 0) go to 7000
C
C Compute transformation matrix from
C Principal to Standard Basis
C
do i = 1,3
C
k = 1 + mod(i,3)
C
do j = 1,3
C
l = 1 + mod(j,3)
C
q(i,j) = qen(i,j)*qen(i,j)
C
q(i+3,j) = qen(i,j)*qen(k,j)
C
q(j,i+3) = qen(j,i)*qen(j,k)*two
C
q(j+3,i+3) = qen(j,i)*qen(l,k) + qen(j,k)*qen(l,i)
C
end do ! j
C
end do ! i
C
C Compute Deviatoric Strains, deviatoric Kirchhoff
C Stresses and Material Moduli in Principal Basis
C
call wlog3d (mu, bpr, epsd, taup, aap, EngW, ipswb3, iwr, ierrl)
if (ierrl .lt. 0) go to 7000
C
C Transform deviatoric Kirchhoff Stresses and Material Moduli
C to Cauchy stresses and spatial material moduli
C
call amat3d (taup, aap, q, Sig, Cst, ipswb3, iwr, ierrl)
if (ierrl .lt. 0) go to 7000
C
C Transform deviatoric Kirchhoff Stresses and Material Moduli
C to Cauchy stresses and spatial material moduli
C
detFi = one / detF
Sig = Sig * detFi
Cst = Cst * detFi
C
C Compute volumetric stresses (pressure) and material moduli
C
iVF = 1
C
call vfnh3d (iVF, detF, detFi, lambda, EngU, Press, Cstv,
& ipswb3, iwr, ierrl)
if (ierrl .lt. 0) go to 7000
C
C Accumulate deviatoric and volumetric contribution to
C spatial stresses and material moduli
C
Sig(1:3) = Sig(1:3) + Press
C
Cstv = Cstv * detF
C
Cst(1,1) = Cst(1,1) - Press + Cstv
Cst(1,2) = Cst(1,2) + Press + Cstv
Cst(1,3) = Cst(1,3) + Press + Cstv
Cst(2,1) = Cst(2,1) + Press + Cstv
Cst(2,2) = Cst(2,2) - Press + Cstv
Cst(2,3) = Cst(2,3) + Press + Cstv
Cst(3,1) = Cst(3,1) + Press + Cstv
Cst(3,2) = Cst(3,2) + Press + Cstv
Cst(3,3) = Cst(3,3) - Press + Cstv
Cst(4,4) = Cst(4,4) - Press
Cst(5,5) = Cst(5,5) - Press
Cst(6,6) = Cst(6,6) - Press
C
C Compute strain energy desnity
C
Engy = EngU + EngW
go to 8000
C
C Error section
C
7000 continue
ierr =-1
go to 7900
7100 continue
ierr =-2
go to 7900
C ----------
7900 continue
write(iwr,9020) '*** ERROR IN SUBROUTINE IEPS3D *** IERR = ', ierr
write(iwr,9010) 'WITH INPUT ARGUMENTS'
write(iwr,9030) 'detF =', detF
write(iwr,9030) 'lambda =', lambda
write(iwr,9030) 'mu =', mu
call rprin0(F , 3, 3, 'F ', iwr)
C
write(iwr,9010) 'WITH OUTPUT ARGUMENTS'
write(iwr,9030) 'Engy =', Engy
call rprin0(Sig, 1, 6,'Sig ', iwr)
call rprin0(Cst, 6, 6,'Cst ', iwr)
go to 8010
C
C Closing section
C
8000 continue
C
if (ipswb3 .gt. 0) then
write(iwr,9010) 'LEAVING SUBROUTINE IEPS3D'
if (ipswb3 .gt. 3) then
write(iwr,9010) 'WITH OUTPUT ARGUMENTS'
write(iwr,9030) 'Engy =', Engy
call rprin0(Sig, 1, 6,'Sig ', iwr)
call rprin0(Cst, 6, 6,'Cst ', iwr)
endif
endif
C
8010 continue
C
return
C
9010 format(3X,A)
9020 format(3X,A,I6)
9030 format(3X,A,E15.6)
C
end

View File

@@ -0,0 +1,183 @@
subroutine liel3d (detF, F, Bmod, Smod, Engy, Sig, Cst,
& ipswb3, iwr, ierr)
C
C ---------------------------------------------------------------------
C
C F E N R I S ROUTINE: LIEL3D NTNU/SINTEF/DNV
C
C PURPOSE:
C LIEL3D: Compute Cauchy stresses and spatial constitutive tensor
C for linear elastic isotropic 3D materials
C
C METHOD:
C Material and spatial stress components are stored as:
C
C | Stress(1) Stress(4) Stress(6) |
C Stress(3,3) = | Stress(4) Stress(2) Stress(5) |
C | Stress(6) Stress(5) Stress(3) |
C
C
C ARGUMENTS INPUT:
C detF - Determinant of the deformation gradient
C F - Deformation gradients
C Bmod - Bulk modulus
C Smod - Shear modulus
C
C ARGUMENTS OUTPUT:
C Engy - Strain energy density
C Sig(6) - Cauchy stresses (spatial stresses)
C Cst(6,6) - Spatial constitutive tensor
C ierr - Error flag
C = 1 : Illegal number of element nodes
C
C COMMON BLOCKS:
C const
C
C PERIPHERAL UNITS:
C None
C
C INTERNAL VARIABLES:
C GLe(6) - Green-Lagrangse strains
C PK2(6) - 2nd Piola-Kirchhoff stresses (material stresses)
C Cmt(6,6) - Material constitutive tensor
C
C PRINT SWITCH:
C ipswb3 = 0 Gives no print
C ipswb3 = 2 Gives enter and leave
C ipswb3 = 3 Gives in addition parameters on input
C ipswb3 = 5 Gives in addition parameters on output
C
C LIMITS:
C
C CODING:
C STANDART ADHERED TO: Computas report no 78-949
C AUTHOR: Kjell Magne Mathisen
C NTNU, Department of structural engineering
C September 14, 2010
C
C INSPECTED:
C
C REVISIONS:
C xx-OCT-2010 K.M.Mathisen NTNU
C Changed XG1-XG8 to XG and UG1-UG8 to UG.
C
C ---------------------------------------------------------------------
C
implicit none
C
integer ipswb3, iwr, ierr
integer i, j, i3
C
real*8 detF, Engy, F(9), Bmod, Smod, Sig(6), Cst(6,6)
real*8 c1, c2
real*8 GLe(6), PK2(6), Cmt(6,6)
C
include 'include/feninc/const.h'
C
C Entry section
C
ierr = 0
C
if (ipswb3 .gt. 0) then
write(iwr,9010) 'ENTERING SUBROUTINE LIEL3D'
if (ipswb3 .gt. 2) then
write(iwr,9010) 'WITH INPUT ARGUMENTS'
write(iwr,9030) 'detF =', detF
write(iwr,9030) 'Bmod =', Bmod
write(iwr,9030) 'Smod =', Smod
call rprin0(F , 3, 3, 'F ', iwr)
endif
endif
C
Cmt = zero
C
C Express the material constitutive tensor
C in terms of the Shear and Bulk modulus
C
c1 = Bmod + four3 * Smod
c2 = Bmod - two3 * Smod
C
do i = 1, 3
j = mod(i,3) + 1
i3 = i + 3
Cmt(i ,i ) = c1
Cmt(i ,j ) = c2
Cmt(j ,i ) = c2
Cmt(i3,i3) = Smod
end do ! i
C
C Compute Green-Lagrange strains from the deformation gradient
C
C E = (F^T * F - 1) / 2
C
GLe(1) = one2 * (F(1)*F(1) + F(2)*F(2) + F(3)*F(3) - one)
GLe(2) = one2 * (F(4)*F(4) + F(5)*F(5) + F(6)*F(6) - one)
GLe(3) = one2 * (F(7)*F(7) + F(8)*F(8) + F(9)*F(9) - one)
GLe(4) = (F(1)*F(4) + F(2)*F(5) + F(3)*F(6))
GLe(5) = (F(4)*F(7) + F(5)*F(8) + F(6)*F(9))
GLe(6) = (F(1)*F(7) + F(2)*F(8) + F(3)*F(9))
C
C Compute 2nd Piola Kirchhoff stresses from product of
C material constitutive tensor and Green-Lagrange strains
C
C PK2 = C_mt * GLe
C
PK2 = matmul(Cmt,GLe)
C
C Compute strain energy density
C
Engy = one2 * dot_product(PK2,GLe)
C
C Push stresses and constitutive tensor forward from
C material to spatial coordinates in current configuration
C
call push3D (detF, F, PK2, Cmt, Sig, Cst, ipswb3, iwr)
C
go to 8000
C
C Error section
C
7000 continue
ierr =-1
go to 7900
7100 continue
ierr =-2
go to 7900
C ----------
7900 continue
write(iwr,9020) '*** ERROR IN SUBROUTINE LIEL3D *** IERR = ', ierr
write(iwr,9010) 'WITH INPUT ARGUMENTS'
write(iwr,9030) 'detF =', detF
write(iwr,9030) 'Bmod =', Bmod
write(iwr,9030) 'Smod =', Smod
call rprin0(F , 3, 3, 'F ', iwr)
C
write(iwr,9010) 'WITH OUTPUT ARGUMENTS'
write(iwr,9030) 'Engy =', Engy
call rprin0(Sig, 1, 6,'Sig ', iwr)
call rprin0(Cst, 6, 6,'Cst ', iwr)
go to 8010
C
C Closing section
C
8000 continue
C
if (ipswb3 .gt. 0) then
write(iwr,9010) 'LEAVING SUBROUTINE LIEL3D'
if (ipswb3 .gt. 3) then
write(iwr,9010) 'WITH OUTPUT ARGUMENTS'
write(iwr,9030) 'Engy =', Engy
call rprin0(Sig, 1, 6,'Sig ', iwr)
call rprin0(Cst, 6, 6,'Cst ', iwr)
endif
endif
C
8010 continue
C
return
C
9010 format(3X,A)
9020 format(3X,A,I6)
9030 format(3X,A,E15.6)
C
end

View File

@@ -0,0 +1,132 @@
subroutine mdma3d (P_bar, P_mix, Sig, DD, ipswb3, iwr)
C
C ---------------------------------------------------------------------
C
C F E N R I S ROUTINE: MDMA3D NTNU/SINTEF/DNV
C
C PURPOSE:
C MDMA3D: Compute finite deformation mixed D arrays
C
C METHOD:
C | D_11 D_12 |
C DD = | |
C | D_21 D_22 |
C where:
C D_11 = 6 x 6 Deviatoric part of matrix
C D_12 = 6 x 1 Coupling part of matrix
C D_21 = 1 x 6 Coupling part of matrix
C D_22 = 1 x 1 Volumetric part of matrix
C
C ARGUMENTS INPUT:
C P_bar - Constitutive pressure
C P_mix - Mixed pressure
C Sig(6) - Constitutive Cauchy stresses
C DD(7,7) - Modified constitutive array
C
C ARGUMENTS OUTPUT:
C DD(7,7) - Mixed material tangent terms
C
C COMMON BLOCKS:
C None
C
C PERIPHERAL UNITS:
C None
C
C INTERNAL VARIABLES:
C None
C
C PRINT SWITCH:
C ipswb3 = 0 Gives no print
C ipswb3 = 2 Gives enter and leave
C ipswb3 = 3 Gives in addition parameters on input
C ipswb3 = 5 Gives in addition parameters on output
C
C LIMITS:
C
C CODING:
C STANDART ADHERED TO: Computas report no 78-949
C AUTHOR: Kjell Magne Mathisen
C NTNU, Department of structural engineering
C September 14, 2010
C
C INSPECTED:
C
C REVISIONS:
C xx-OCT-2010 K.M.Mathisen NTNU
C Changed XG1-XG8 to XG and UG1-UG8 to UG.
C
C ---------------------------------------------------------------------
C
implicit none
C
integer ipswb3, iwr, j
real*8 DD(7,7), P_bar, P_mix, Sig(6)
real*8 fac1, Sig_D(6)
C
include 'include/feninc/const.h'
C
C Entry section
C
if (ipswb3 .gt. 0) then
write(iwr,9010) 'ENTERING SUBROUTINE MDMA3D'
if (ipswb3 .gt. 2) then
write(iwr,9010) 'WITH INPUT ARGUMENTS'
write(iwr,9030) 'P_bar =', P_bar
write(iwr,9030) 'P_mix =', P_mix
call rprin0(Sig, 1, 6, 'Sig ', iwr)
call rprin0(DD , 7, 7, 'DD ', iwr)
endif
endif
C
C
C Compute deviatoric stress
C
Sig_D(1:3) = two3 *(Sig(1:3) - P_bar)
Sig_D(4:6) = two3 * Sig(4:6)
C
C D_11: B_matrix part
C
fac1 = P_mix - two3 * P_bar
C
DD(1:3,1:3) = DD(1:3,1:3) + fac1
C
do j = 1,6
DD(1:3,j) = DD(1:3,j) - Sig_D(j)
DD(j,1:3) = DD(j,1:3) - Sig_D(j)
end do
C
fac1 = P_bar - P_mix
C
do j = 1,3
DD(j ,j ) = DD(j ,j ) + fac1 * two
DD(j+3,j+3) = DD(j+3,j+3) + fac1
end do ! j
C
C D_12: Coupling matrix with
C
DD(7,1:6) = DD(7,1:6) + Sig_D
DD(1:6,7) = DD(1:6,7) + Sig_D
C
C D_22: Volumetric part
C
DD(7,7) = DD(7,7) - one3 * P_bar
C
C
C Closing section
C
8000 continue
if (ipswb3 .gt. 0) then
write(iwr,9010) 'LEAVING SUBROUTINE MDMA3D'
if (ipswb3 .gt. 3) then
write(iwr,9010) 'WITH OUTPUT ARGUMENTS'
call rprin0(DD , 7, 7, 'DD ', iwr)
endif
call flush(iwr)
endif
C
return
C
9010 format(3X,A)
9030 format(3X,A,E15.6)
C
end

View File

@@ -0,0 +1,274 @@
subroutine monh3d (iVF, detF, F, Bmod, Smod, Engy, Sig, Cst,
& ipswb3, iwr, ierr)
C
C ---------------------------------------------------------------------
C
C F E N R I S ROUTINE: MONH3D NTNU/SINTEF/DNV
C
C PURPOSE:
C MONH3D: Compute Cauchy stresses and spatial constitutive tensor
C for a modified NeoHookean isotropic 3D material
C
C METHOD:
C Energy function: Compressible Neohookean model
C with J_2/3 regularization
C _ __ _ _ __
C W(J,be) = U(J) + 0.5*mu*(J^(-2/3)be:1 - 3)
C
C Volumetric function type: U_1 = K*[ 0.25*( J^2 - 1 )-0.5*ln(J)]
C U_2 = K*[ 0.5*( J - 1 )^2 ]
C U_3 = K*[ 0.5*( ln(J) )^2 ]
C U_4 = K*[ 2.0*(J - 1 - ln(J)) ]
C Free energy function:
C
C Psi = U(J) -3*alpha*J*U'(J)*[T - T_ref]
C
C up = partial_J ( Psi )
C upp = [ partial_J ( J partial_J Psi ) - up ]/J
C = partial^2_J ( Psi )/JC
C
C
C ARGUMENTS INPUT:
C iVF - Volumetric function type
C detF - Determinant of the deformation gradient
C F - Deformation gradients
C Bmod - Bulk modulus
C Smod - Shear modulus
C
C ARGUMENTS OUTPUT:
C Engy - Strain energy density
C Sig(6) - Cauchy stresses (spatial stresses)
C Cst(6,6) - Spatial constitutive tensor
C ierr - Error flag
C = 1 : Illegal number of element nodes
C
C COMMON BLOCKS:
C const
C
C PERIPHERAL UNITS:
C None
C
C INTERNAL VARIABLES:
C nSTR - Number of stress/strain components
C bm(6) - Modified Left Cauchy-Green deformation tensor
C
C
C PRINT SWITCH:
C ipswb3 = 0 Gives no print
C ipswb3 = 2 Gives enter and leave
C ipswb3 = 3 Gives in addition parameters on input
C ipswb3 = 5 Gives in addition parameters on output
C
C LIMITS:
C
C CODING:
C STANDART ADHERED TO: Computas report no 78-949
C AUTHOR: Kjell Magne Mathisen
C NTNU, Department of structural engineering
C September 14, 2010
C
C INSPECTED:
C
C REVISIONS:
C xx-OCT-2010 K.M.Mathisen NTNU
C Changed XG1-XG8 to XG and UG1-UG8 to UG.
C
C ---------------------------------------------------------------------
C
implicit none
C
integer iVF, ipswb3, iwr, ierr
integer i, j, ierrl
C
real*8 detF, Engy, F(3,3), Bmod, Smod, Sig(6), Cst(6,6)
real*8 c1, c2, c3, Cstv, bm(6), detFi, Press, trbm3
C
include 'include/feninc/const.h'
C
C Entry section
C
ierr = 0
C
if (ipswb3 .gt. 0) then
write(iwr,9010) 'ENTERING SUBROUTINE MONH3D'
if (ipswb3 .gt. 2) then
write(iwr,9010) 'WITH INPUT ARGUMENTS'
write(iwr,9020) 'iVF =', iVF
write(iwr,9030) 'detF =', detF
write(iwr,9030) 'Bmod =', Bmod
write(iwr,9030) 'Smod =', Smod
call rprin0(F , 3, 3, 'F ', iwr)
endif
endif
C
C Initialize spatial constitutive tensor
C
Cst = zero
C
C Compute Left Cauchy-Green deformation tensor
C
C b = F * F^T
C
bm(1) = F(1,1)*F(1,1) + F(1,2)*F(1,2) + F(1,3)*F(1,3)
bm(2) = F(2,1)*F(2,1) + F(2,2)*F(2,2) + F(2,3)*F(2,3)
bm(3) = F(3,1)*F(3,1) + F(3,2)*F(3,2) + F(3,3)*F(3,3)
bm(4) = F(1,1)*F(2,1) + F(1,2)*F(2,2) + F(1,3)*F(2,3)
bm(5) = F(2,1)*F(3,1) + F(2,2)*F(3,2) + F(2,3)*F(3,3)
bm(6) = F(1,1)*F(3,1) + F(1,2)*F(3,2) + F(1,3)*F(3,3)
C
if (ipswb3 .gt. 3) then
call rprin0(bm, 1, 6,'bm ', iwr)
endif
C
C Compute modified Left Cauchy-Green deformation tensor
C _
C bm = b = J^(-2/3) * b
C
detFi = one / detF
bm = bm * (detFi ** two3)
C
if (ipswb3 .gt. 3) then
call rprin0(bm, 1, 6,'bm ', iwr)
endif
C
C Compute one third of the trace of the modified
C Left Cauchy-Green deformation tensor
C _
C trbm3 = tr( b ) / 3
C
trbm3 = (bm(1) + bm(2) + bm(3)) * one3
bm(1) = bm(1) - trbm3
bm(2) = bm(2) - trbm3
bm(3) = bm(3) - trbm3
C
if (ipswb3 .gt. 3) then
call rprin0(bm, 1, 6,'bm ', iwr)
endif
C
C Compute deviatoric part of the constitutive tensor
C
C Part 1: _ _ _
C Rank one update: -2/3 G * ( b x g + g x b ) / J
C
c1 = two3 * Smod
C
do i = 1, 6
c2 = c1 * bm(i)
do j = 1,3
Cst(i,j) = Cst(i,j) - c2
Cst(j,i) = Cst(j,i) - c2
end do ! j
end do ! i
C
C Part 2: __ _
C Deviatoric term 2 mu [ I - 1/3 g x g ] / J
C
c1 = Smod * trbm3
c2 = c1 + c1
c3 = c2 * one3
C
do i = 1,3
Cst(i ,i ) = Cst(i ,i ) + c2
Cst(i+3,i+3) = Cst(i+3,i+3) + c1
do j = 1,3
Cst(i ,j ) = Cst(i ,j ) - c3
end do ! j
end do ! i
C
C Compute deviatoric Kirchhoff stress tensor
C _ _
C tau_dev = b - tr( b ) / 3
C
Sig = Smod * bm
C
C Compute spatial deviatoric (isochoric)
C stresses and material moduli
C
Sig = Sig * detFi
Cst = Cst * detFi
C
if (ipswb3 .gt. 3) then
write(iwr,9030) 'Press =', Press
call rprin0(Sig, 1, 6,'Sig ', iwr)
call rprin0(Cst, 6, 6,'Cst ', iwr)
endif
C
C Compute volumetric stresses (presssure) and material moduli
C according to the volumetric function type
C
call vfnh3d (iVF, detF, detFi, Bmod, Engy, Press, Cstv,
& ipswb3, iwr, ierrl)
if (ierrl .lt. 0) go to 7000
C
Sig(1:3) = Sig(1:3) + Press
C
Cstv = Cstv * detF
C
Cst(1,1) = Cst(1,1) - Press + Cstv
Cst(1,2) = Cst(1,2) + Press + Cstv
Cst(1,3) = Cst(1,3) + Press + Cstv
Cst(2,1) = Cst(2,1) + Press + Cstv
Cst(2,2) = Cst(2,2) - Press + Cstv
Cst(2,3) = Cst(2,3) + Press + Cstv
Cst(3,1) = Cst(3,1) + Press + Cstv
Cst(3,2) = Cst(3,2) + Press + Cstv
Cst(3,3) = Cst(3,3) - Press + Cstv
Cst(4,4) = Cst(4,4) - Press
Cst(5,5) = Cst(5,5) - Press
Cst(6,6) = Cst(6,6) - Press
C
C Compute strain energy density
C
Engy = Engy + onep5 * Smod * (trbm3 - one)
go to 8000
C
C Error section
C
7000 continue
ierr =-1
go to 7900
7100 continue
ierr =-2
go to 7900
C ----------
7900 continue
write(iwr,9020) '*** ERROR IN SUBROUTINE MONH3D *** IERR = ', ierr
write(iwr,9010) 'WITH INPUT ARGUMENTS'
write(iwr,9020) 'iVF =', iVF
write(iwr,9030) 'detF =', detF
write(iwr,9030) 'Bmod =', Bmod
write(iwr,9030) 'Smod =', Smod
call rprin0(F , 3, 3, 'F ', iwr)
C
write(iwr,9010) 'WITH OUTPUT ARGUMENTS'
write(iwr,9030) 'Engy =', Engy
call rprin0(Sig, 1, 6,'Sig ', iwr)
call rprin0(Cst, 6, 6,'Cst ', iwr)
go to 8010
C
C Closing section
C
8000 continue
C
if (ipswb3 .gt. 0) then
write(iwr,9010) 'LEAVING SUBROUTINE MONH3D'
if (ipswb3 .gt. 3) then
write(iwr,9010) 'WITH OUTPUT ARGUMENTS'
write(iwr,9030) 'Engy =', Engy
call rprin0(Sig, 1, 6,'Sig ', iwr)
call rprin0(Cst, 6, 6,'Cst ', iwr)
endif
endif
C
8010 continue
C
return
C
9010 format(3X,A)
9020 format(3X,A,I6)
9030 format(3X,A,E15.6)
C
end

View File

@@ -0,0 +1,121 @@
subroutine pcst3d (aa,dd,ipswb3,iwr)
C
C ---------------------------------------------------------------------
C
C F E N R I S ROUTINE: PCST3D NTNU/SINTEF/DNV
C
C PURPOSE:
C PCST3D: Project 6x6 constitutive tensor (AA) onto a 7x7 equivalent
C constitutive tensor for mixed model (DD)
C
C METHOD:
C | D_11 D_12 |
C DD = | |
C | D_21 D_22 |
C where:
C D_11 = 6 x 6 Deviatoric part of matrix
C D_12 = 6 x 1 Coupling part of matrix
C D_21 = 1 x 6 Coupling part of matrix
C D_22 = 1 x 1 Volumetric part of matrix
C
C
C ARGUMENTS INPUT:
C aa(6,6) - Material tangent matrix (based on F)
C
C ARGUMENTS OUTPUT:
C dd(7,7) - Mixed material tangent for stiffness computations
C
C COMMON BLOCKS:
C None
C
C PERIPHERAL UNITS:
C None
C
C INTERNAL VARIABLES:
C None
C
C PRINT SWITCH:
C ipswb3 = 0 Gives no print
C ipswb3 = 2 Gives enter and leave
C ipswb3 = 3 Gives in addition parameters on input
C ipswb3 = 5 Gives in addition parameters on output
C
C LIMITS:
C
C CODING:
C STANDART ADHERED TO: Computas report no 78-949
C AUTHOR: Kjell Magne Mathisen
C NTNU, Department of structural engineering
C September 14, 2010
C
C INSPECTED:
C
C REVISIONS:
C xx-OCT-2010 K.M.Mathisen NTNU
C Changed XG1-XG8 to XG and UG1-UG8 to UG.
C
C ---------------------------------------------------------------------
C
implicit none
C
integer ipswb3, iwr, j
real*8 aa(6,6), dd(7,7)
C
include 'include/feninc/const.h'
C
C
C Entry section
C
if (ipswb3 .gt. 0) then
write(iwr,9010) 'ENTERING SUBROUTINE PCST3D'
if (ipswb3 .gt. 2) then
write(iwr,9010) 'WITH INPUT ARGUMENTS'
call rprin0(aa, 6, 6, 'aa ', iwr)
endif
endif
C
C
C Load moduli from constitution
C
dd(1:6,1:6) = aa
C
C Compute left and right multiples with trace
C
dd(1:6,7) = (aa(:,1) + aa(:,2) + aa(:,3))*one3
dd(7,1:6) = (aa(1,:) + aa(2,:) + aa(3,:))*one3
C
C Convert upper 6 x 6 to a deviatoric D_11
C
do j = 1, 3
dd(1:6,j) = dd(1:6,j) - dd(1:6,7)
dd(j,1:6) = dd(j,1:6) - dd(7,1:6)
end do
C
C Form last term, D_22
C
dd(7,7) = (dd(1,7) + dd(2,7) + dd(3,7))*one3
C
C Final update to form D_12 and D_21
C
dd(1:3, 7) = dd(1:3, 7) - dd(7,7)
dd( 7,1:3) = dd( 7,1:3) - dd(7,7)
dd(1:3,1:3) = dd(1:3,1:3) + dd(7,7)
C
C
C Closing section
C
8000 continue
if (ipswb3 .gt. 0) then
write(iwr,9010) 'LEAVING SUBROUTINE PCST3D'
if (ipswb3 .gt. 3) then
write(iwr,9010) 'WITH OUTPUT ARGUMENTS'
call rprin0(dd, 7, 7, 'dd ', iwr)
endif
call flush(iwr)
endif
C
return
C
9010 format(3X,A)
C
end

View File

@@ -0,0 +1,162 @@
subroutine push3d (detF, F, PK2, Cmt, Sig, Cst, ipswb3, iwr)
C
C ---------------------------------------------------------------------
C
C F E N R I S ROUTINE: PUSH3D NTNU/SINTEF/DNV
C
C PURPOSE:
C PUSH3D: Push stresses and constitutive tensor forward from
C material to spatial coordinates in current configuration
C
C METHOD:
C Material and spatial stress components are stored as:
C
C | Stress(1) Stress(4) Stress(6) |
C Stress(3,3) = | Stress(4) Stress(2) Stress(5) |
C | Stress(6) Stress(5) Stress(3) |
C
C ARGUMENTS INPUT:
C detF - Determinant of deformation gradient
C F(3,3) - Deformation gradient
C PK2(6) - 2nd Piola-Kirchhoff stresses (material stresses)
C Cmt(6,6) - Material constitutive tensor
C
C ARGUMENTS OUTPUT:
C Sig(6) - Cauchy stresses (spatial stresses)
C Cst(6,6) - Spatial constitutive tensor
C
C COMMON BLOCKS:
C const
C
C PERIPHERAL UNITS:
C None
C
C INTERNAL VARIABLES:
C Tmat(6,6) - Transformation matrix
C
C PRINT SWITCH:
C ipswb3 = 0 Gives no print
C ipswb3 = 2 Gives enter and leave
C ipswb3 = 3 Gives in addition parameters on input
C ipswb3 = 5 Gives in addition parameters on output
C
C LIMITS:
C
C CODING:
C STANDART ADHERED TO: Computas report no 78-949
C AUTHOR: Kjell Magne Mathisen
C NTNU, Department of structural engineering
C September 14, 2010
C
C INSPECTED:
C
C REVISIONS:
C xx-OCT-2010 K.M.Mathisen NTNU
C Changed XG1-XG8 to XG and UG1-UG8 to UG.
C
C ---------------------------------------------------------------------
C
implicit none
C
integer ipswb3, iwr
integer i, j, i1(6), i2(6)
real*8 detF, F(3,3), PK2(6), Cmt(6,6), Sig(6), Cst(6,6)
real*8 detFi, Tmat(6,6)
C
include 'include/feninc/const.h'
C
data i1 /1,2,3,1,2,3/
data i2 /1,2,3,2,3,1/
C
C Entry section
C
if (ipswb3 .gt. 0) then
write(iwr,9010) 'ENTERING SUBROUTINE PUSH3D'
if (ipswb3 .gt. 2) then
write(iwr,9010) 'WITH INPUT ARGUMENTS'
write(iwr,9020) 'detF =', detF
call rprin0(F , 3, 3,'F ',iwr)
call rprin0(PK2, 1, 6,'PK2 ',iwr)
call rprin0(Cmt, 6, 6,'Cmt ',iwr)
endif
endif
C
C Form transformation array for a 4th rank tensor in matrix form
C
C Tmat(a,b) = F(i,I)*F(j,J) : a -> I,J ; b -> i,j
C a,b | 1 2 3 4 5 6
C ------+-----------------------------
C (I,J) | 1,1 2,2 3,3 1,2 2,3 3,1
C or (i,j) | 2,1 3,2 1,3
C
do i = 1,3
do j = 1,3
Tmat(i,j) = F(i1(j),i1(i))*F(i2(j),i2(i))
end do ! j
do j = 4,6
Tmat(i,j) = (F(i1(j),i1(i))*F(i2(j),i2(i))
& + F(i2(j),i2(i))*F(i1(j),i1(i)))*one2
end do ! j
end do ! i
do i = 4,6
do j = 1,3
Tmat(i,j) = F(i1(j),i1(i))*F(i2(j),i2(i))
& + F(i2(j),i2(i))*F(i1(j),i1(i))
end do ! j
do j = 4,6
Tmat(i,j) = (F(i1(j),i1(i))*F(i2(j),i2(i))
& + F(i2(j),i1(i))*F(i1(j),i2(i))
& + F(i1(j),i2(i))*F(i2(j),i1(i))
& + F(i2(j),i2(i))*F(i1(j),i1(i)))*one2
end do ! j
end do ! i
C
C
C Compute matrix product: Cst = Tmat^t * Cmt * Tmat
C
Cst = matmul(transpose(Tmat),matmul(Cmt,Tmat))
C
C Push forward material stress (PK2) to current configuration
C
C Sig(i,j) = F(i,k)*PK2(k,l)*F(j,l)/detf
C
C
C Tmat = F^T * PK2
C
do i = 1,3
Tmat(i,1) = F(i,1)*PK2(1) + F(i,2)*PK2(4) + F(i,3)*PK2(6)
Tmat(i,2) = F(i,1)*PK2(4) + F(i,2)*PK2(2) + F(i,3)*PK2(5)
Tmat(i,3) = F(i,1)*PK2(6) + F(i,2)*PK2(5) + F(i,3)*PK2(3)
end do ! i
C
C Sig = Tmat*F
C
Sig(1) = Tmat(1,1)*F(1,1) + Tmat(1,2)*F(1,2) + Tmat(1,3)*F(1,3)
Sig(2) = Tmat(2,1)*F(2,1) + Tmat(2,2)*F(2,2) + Tmat(2,3)*F(2,3)
Sig(3) = Tmat(3,1)*F(3,1) + Tmat(3,2)*F(3,2) + Tmat(3,3)*F(3,3)
Sig(4) = Tmat(1,1)*F(2,1) + Tmat(1,2)*F(2,2) + Tmat(1,3)*F(2,3)
Sig(5) = Tmat(2,1)*F(3,1) + Tmat(2,2)*F(3,2) + Tmat(2,3)*F(3,3)
Sig(6) = Tmat(3,1)*F(1,1) + Tmat(3,2)*F(1,2) + Tmat(3,3)*F(1,3)
C
detFi = one / detf
Cst = detFi * Cst
Sig = detFi * Sig
C
C Closing section
C
if (ipswb3 .gt. 0) then
write(iwr,9010) 'LEAVING SUBROUTINE PUSH3D'
if (ipswb3 .gt. 3) then
write(iwr,9010) 'WITH OUTPUT ARGUMENTS'
call rprin0(Sig, 1, 6,'Sig ',iwr)
call rprin0(Cst, 6, 6,'Cst ',iwr)
endif
endif
C
return
C
9010 format(3X,A)
9020 format(3X,A,E15.6)
C
end

View File

@@ -0,0 +1,93 @@
SUBROUTINE RPRIN0 (A,MA,NA,CHHEAD,IW)
INTEGER MA,NA, IW
DOUBLE PRECISION A(MA,* )
CHARACTER CHHEAD*(*)
C
C @(#)rprin0.f 1.3 90/11/19 Copyright (c) 1989 by A.S Veritas Research.
C-----------------------------------------------------------------------
C
C F E N R I S ROUTINE : RPRIN0 NTH/SINTEF/DNV
C
C PURPOSE:
C RPRIN0: Print the real matrix A(MA,NA).
C
C METHOD:
C THE MATRIX IS PRINTED AS A TWO-DIMENSIONAL MATRIX.ROWS AND COLUMNS
C ARE NUMBERED,AND EACH ELEMENT IS PRINTED WITH THE FORMAT 1PE14.6 ,
C I.E. 7 SIGNIFICANTS DIGITS.
C A HEADING ( E.G. THE SYMBOL NAME "A" OF THE MATRIX ) CONSISTING OF
C UP TO 6 CHARACTERS IS PRINTED.
C EXAMPLE OF A CALL TO RPRIN0: CALL RPRIN0(A,MA,NA,'A ',IW)
C
C ARGUMENTS INPUT:
C A - MATRIX TO BE PRINTED
C MA - NUMBER OF ROWS
C NA - NUMBER OF COLUMNS
C CHHEAD - HEADING, E.G. "NAME" OF MATRIX
C IW - LINE PRINTER UNIT NUMBER
C
C ARGUMENTS OUTPUT:
C NONE
C
C COMMON BLOCKS:
C NONE
C
C PERIPHERAL UNITS:
C NONE
C
C LIMITATIONS:
C NONE
C
C CODING:
C STANDARD ADHERED TO: COMPUTAS REPORT NO 78-949
C CODED: DECEMBER 5, 1974 K.BELL SINTEF - DIV. OF STRUCT. ENG.
C
C REVISIONS:
C DATE: AUGUST 27,1980 ALF ENGSETH SINTEF - DIV. OF STRUCT. ENG.
C DATE: DECEMBER 20,1982 MAGNE NYGAARD VERITAS
C 12-MAR-87 M.K.Nygaard VERITEC 86.1-0
C Changed format 9100
C 2-DEC-87 K.M.Mathisen NTH 5.2-04
C Included test on the value of NA and changed dimensioning of A.
C 16-NOV-88 Ole Stroem / Robert O. Bjaerum, Veritas Research
C Any dimension of CHHEAD OK. Please keep it under 120.
C
C-----------------------------------------------------------------------
C
INTEGER NCOL,J,I,K,L
PARAMETER ( NCOL = 9 )
C
C EXECUTABLE SECTION
C
C PROGRAM LOGIC SECTION
C
IF (NA .LE. 0 .OR. MA .LE. 0) GO TO 8000
C
C PRINT HEADING
C
WRITE(IW,9050) CHHEAD
C
C PRINT MATRIX
C
DO 4000 J=1,NA,NCOL
K = J+(NCOL-1)
IF ( K .GT. NA ) THEN
K = NA
END IF
C
WRITE(IW,9000) (L,L=J,K)
C
DO 3000 I=1,MA
WRITE(IW,9100) I,(A(I,L),L=J,K)
3000 CONTINUE
4000 CONTINUE
C
8000 CONTINUE
C
RETURN
C
9000 FORMAT ( / 3X , 9( I9 , 5X ) )
9050 FORMAT ( / 5X , A , ' = ' )
9100 FORMAT ( 1X , I3 , 1P , 9E14.6 )
C
END

View File

@@ -0,0 +1,210 @@
subroutine stnh3d (iVF, detF, F, lambda, mu, Engy, Sig, Cst,
& ipswb3, iwr, ierr)
C
C ---------------------------------------------------------------------
C
C F E N R I S ROUTINE: STNH3D NTNU/SINTEF/DNV
C
C PURPOSE:
C MONH3D: Compute Cauchy stresses and spatial constitutive tensor
C for a standard NeoHookean isotropic 3D material
C
C METHOD:
C Energy function: Standard Neohookean model
C
C W(J,b) = lambda*U(J) + 0.5*mu*(I_c - 3 - 2*ln(J))
C
C Sig = [mu*(b - 1) + lambda*Press]/detF
C = mu*(b - 1)/detF + lambda*partial_U/partial_J
C
C Volumetric function type: U_1 = K*[ 0.25*( J^2 - 1 )-0.5*ln(J)]
C U_2 = K*[ 0.5*( J - 1 )^2 ]
C U_3 = K*[ 0.5*( ln(J) )^2 ]
C U_4 = K*[ 2.0*(J - 1 - ln(J)) ]
C Free energy function:
C
C Psi = U(J) -3*alpha*J*U'(J)*[T - T_ref]
C
C up = partial_J ( Psi )
C upp = [ partial_J ( J partial_J Psi ) - up ]/J
C = partial^2_J ( Psi )/JC
C
C
C ARGUMENTS INPUT:
C iVF - Volumetric function type
C detF - Determinant of the deformation gradient
C F - Deformation gradients
C lambda - Lame's constant (=K-2*G/3)
C mu - Shear modulus
C
C ARGUMENTS OUTPUT:
C Engy - Strain energy density
C Sig(6) - Cauchy stresses (spatial stresses)
C Cst(6,6) - Spatial constitutive tensor
C ierr - Error flag
C = 1 : Illegal number of element nodes
C
C COMMON BLOCKS:
C const
C
C PERIPHERAL UNITS:
C None
C
C INTERNAL VARIABLES:
C b(6) - Left Cauchy-Green deformation tensor
C
C
C PRINT SWITCH:
C ipswb3 = 0 Gives no print
C ipswb3 = 2 Gives enter and leave
C ipswb3 = 3 Gives in addition parameters on input
C ipswb3 = 5 Gives in addition parameters on output
C
C LIMITS:
C
C CODING:
C STANDART ADHERED TO: Computas report no 78-949
C AUTHOR: Kjell Magne Mathisen
C NTNU, Department of structural engineering
C September 14, 2010
C
C INSPECTED:
C
C REVISIONS:
C xx-OCT-2010 K.M.Mathisen NTNU
C Changed XG1-XG8 to XG and UG1-UG8 to UG.
C
C ---------------------------------------------------------------------
C
implicit none
C
integer iVF, ipswb3, iwr, ierr
integer i, i3, ierrl
C
real*8 detF, Engy, F(3,3), lambda, mu, Sig(6), Cst(6,6)
real*8 c1, c2, b(6), Cstv, detFi, Press
C
include 'include/feninc/const.h'
C
C Entry section
C
ierr = 0
C
if (ipswb3 .gt. 0) then
write(iwr,9010) 'ENTERING SUBROUTINE STNH3D'
if (ipswb3 .gt. 2) then
write(iwr,9010) 'WITH INPUT ARGUMENTS'
write(iwr,9020) 'iVF =', iVF
write(iwr,9030) 'detF =', detF
write(iwr,9030) 'lambda =', lambda
write(iwr,9030) 'mu =', mu
call rprin0(F , 3, 3, 'F ', iwr)
endif
endif
C
C Compute Left Cauchy-Green deformation tensor
C
C b = F * F^T
C
b(1) = F(1,1)*F(1,1) + F(1,2)*F(1,2) + F(1,3)*F(1,3)
b(2) = F(2,1)*F(2,1) + F(2,2)*F(2,2) + F(2,3)*F(2,3)
b(3) = F(3,1)*F(3,1) + F(3,2)*F(3,2) + F(3,3)*F(3,3)
b(4) = F(1,1)*F(2,1) + F(1,2)*F(2,2) + F(1,3)*F(2,3)
b(5) = F(2,1)*F(3,1) + F(2,2)*F(3,2) + F(2,3)*F(3,3)
b(6) = F(1,1)*F(3,1) + F(1,2)*F(3,2) + F(1,3)*F(3,3)
C
C Compute volumetric stresses (pressure) and material moduli
C according to the volumetric function type
C
detFi = one / detF
C
call vfnh3d (iVF, detF, detFi, lambda, Engy, Press, Cstv,
& ipswb3, iwr, ierrl)
if (ierrl .lt. 0) go to 7000
C
C Initialize spatial constitutive tensor
C
Cst = zero
C
C Accumulate deviatoric and volumetric contribution to
C spatial stresses and material moduli
C
c1 = mu * detFi
C c2 = two * mu
c2 = two * c1
Cstv = Cstv * detF
C
C Diagonal terms
C
do i = 1, 3
i3 = i + 3
Sig(i) = c1 * b(i) - c1 + Press
Sig(i3) = c1 * b(i3)
Cst(i ,i ) = c2 - Press + Cstv
Cst(i3,i3) = c1 - Press
end do ! i
C
C Off-diagonal terms
C
c1 = Press + Cstv
C
Cst(1,2) = c1
Cst(1,3) = c1
Cst(2,1) = c1
Cst(2,3) = c1
Cst(3,1) = c1
Cst(3,2) = c1
C
C Compute strain energy desnity
C
Engy = Engy + mu*(one2*(b(1)+b(2)+b(3)) - onep5 - log(abs(detF)))
C
go to 8000
C
C Error section
C
7000 continue
ierr =-1
go to 7900
7100 continue
ierr =-2
go to 7900
C ----------
7900 continue
write(iwr,9020) '*** ERROR IN SUBROUTINE STNH3D *** IERR = ', ierr
write(iwr,9010) 'WITH INPUT ARGUMENTS'
write(iwr,9020) 'iVF =', iVF
write(iwr,9030) 'detF =', detF
write(iwr,9030) 'lambda =', lambda
write(iwr,9030) 'mu =', mu
call rprin0(F , 3, 3, 'F ', iwr)
C
write(iwr,9010) 'WITH OUTPUT ARGUMENTS'
write(iwr,9030) 'Engy =', Engy
call rprin0(Sig, 1, 6,'Sig ', iwr)
call rprin0(Cst, 6, 6,'Cst ', iwr)
go to 8010
C
C Closing section
C
8000 continue
C
if (ipswb3 .gt. 0) then
write(iwr,9010) 'LEAVING SUBROUTINE STNH3D'
if (ipswb3 .gt. 3) then
write(iwr,9010) 'WITH OUTPUT ARGUMENTS'
write(iwr,9030) 'Engy =', Engy
call rprin0(Sig, 1, 6,'Sig ', iwr)
call rprin0(Cst, 6, 6,'Cst ', iwr)
endif
endif
C
8010 continue
C
return
C
9010 format(3X,A)
9020 format(3X,A,I6)
9030 format(3X,A,E15.6)
C
end

View File

@@ -0,0 +1,187 @@
subroutine vfnh3d (iVF, detF, detFi, lambda, U, Up, Upp,
& ipswb3, iwr, ierr)
C
C ---------------------------------------------------------------------
C
C F E N R I S ROUTINE: VFNH3D NTNU/SINTEF/DNV
C
C PURPOSE:
C VFNH3D: Compute volumetric stresses (pressure) and material moduli
C according to the volumetric function type
C
C METHOD:
C Volumetric function type: U_1 = K*[ 0.25*( J^2 - 1 )-0.5*ln(J)]
C U_2 = K*[ 0.5*( J - 1 )^2 ]
C U_3 = K*[ 0.5*( ln(J) )^2 ]
C U_4 = K*[ 2.0*(J - 1 - ln(J)) ]
C Free energy function:
C U = U(J)
C Up = partial_J ( Psi )
C Upp = [ partial_J ( J partial_J Psi ) - up ]/J
C = partial^2_J ( Psi )/JC
C
C
C ARGUMENTS INPUT:
C iVF - Volumetric function type
C detF - Determinant of the deformation gradient
C detFi - Inverse of the determinant of the Jacobian gradient
C lambda - Lame's constant
C
C ARGUMENTS OUTPUT:
C U - Strain energy density
C Up - Mean stress (pressure)
C Upp - Volumetric constitutive tensor
C ierr - Error flag
C = 1 : Illegal number of element nodes
C
C COMMON BLOCKS:
C const
C
C PERIPHERAL UNITS:
C None
C
C INTERNAL VARIABLES:
C U = U(J)
C Up = partial_J ( Psi )
C
C PRINT SWITCH:
C ipswb3 = 0 Gives no print
C ipswb3 = 2 Gives enter and leave
C ipswb3 = 3 Gives in addition parameters on input
C ipswb3 = 5 Gives in addition parameters on output
C
C LIMITS:
C
C CODING:
C STANDART ADHERED TO: Computas report no 78-949
C AUTHOR: Kjell Magne Mathisen
C NTNU, Department of structural engineering
C September 14, 2010
C
C INSPECTED:
C
C REVISIONS:
C xx-OCT-2010 K.M.Mathisen NTNU
C Changed XG1-XG8 to XG and UG1-UG8 to UG.
C
C ---------------------------------------------------------------------
C
implicit none
C
integer iVF, ipswb3, iwr, ierr
C
real*8 detF, detFi, lambda,
& detF2, detFi2, detFm1, U, Up, Upp
C
include 'include/feninc/const.h'
C
C Entry section
C
ierr = 0
C
if (ipswb3 .gt. 0) then
write(iwr,9010) 'ENTERING SUBROUTINE VFNH3D'
if (ipswb3 .gt. 2) then
write(iwr,9010) 'WITH INPUT ARGUMENTS'
write(iwr,9020) 'iVF =', iVF
write(iwr,9030) 'detF =', detF
write(iwr,9030) 'detFi =', detFi
write(iwr,9030) 'lambda =', lambda
endif
endif
C
C Compute volumetric stresses (pressure) and material moduli
C according to the volumetric function type
C
if (iVF . eq. 1) then
C
C Volumetric function type 1:
C U(J) = lambda*0.25*(J^2 - 1 - 2*(log J))
C
detF2 = detF**2
detFi2 = detFi**2
U = one4 * lambda * ( detF2 - one - two*log(abs(detF)) )
Up = one2 * lambda * ( detF - detFi )
Upp = one2 * lambda * ( one + detFi2 )
C
else if (iVF .eq. 2) then
C
C Volumetric function type 2:
C U(J) = lambda*0.5*(J-1)^2
C
detFm1 = detF - one
U = one2 * lambda * (detFm1)**2
Up = lambda * detFm1
Upp = lambda
C
else if (IVF .eq. 3) then
C
C Volumetric function type 3:
C U(J) = lambda*0.5*(log J)^2
C
U = one2 * lambda * log(abs(detF))**2
Up = lambda * log(abs(detF)) * detFi
Upp = ( lambda * detFi - Up ) * detFi
C
else if (iVF .eq. 4) then
C
C Volumetric function type 4:
C U(J) = lambda*2.0*(J - 1 - log J)
C
detFi2 = detFi**2
U = two * lambda * ( detF - one - log(abs(detF)) )
Up = lambda * ( two - detFi )
Upp = lambda * detFi2
C
C Illegal volumetric function
C
else
go to 7000
end if
go to 8000
C
C Error section
C
7000 continue
ierr =-1
go to 7900
7100 continue
ierr =-2
go to 7900
C ----------
7900 continue
write(iwr,9020) '*** ERROR IN SUBROUTINE VFNH3D *** IERR = ', ierr
write(iwr,9010) 'WITH INPUT ARGUMENTS'
write(iwr,9030) 'detF =', detF
write(iwr,9030) 'detFi =', detFi
write(iwr,9030) 'lambda =', lambda
C
write(iwr,9010) 'WITH OUTPUT ARGUMENTS'
write(iwr,9030) 'U =', U
write(iwr,9030) 'Up =', Up
write(iwr,9030) 'Upp =', Upp
go to 8010
C
C Closing section
C
8000 continue
C
if (ipswb3 .gt. 0) then
write(iwr,9010) 'LEAVING SUBROUTINE VFNH3D'
if (ipswb3 .gt. 3) then
write(iwr,9010) 'WITH OUTPUT ARGUMENTS'
write(iwr,9030) 'U =', U
write(iwr,9030) 'Up =', Up
write(iwr,9030) 'Upp =', Upp
endif
endif
C
8010 continue
C
return
C
9010 format(3X,A)
9020 format(3X,A,I6)
9030 format(3X,A,E15.6)
C
end

View File

@@ -0,0 +1,184 @@
subroutine wlog3d (mu, bpr, epsd, taup, aap, w,
& ipswb3, iwr, ierr)
C
C ---------------------------------------------------------------------
C
C F E N R I S ROUTINE: WLOG3D NTNU/SINTEF/DNV
C
C PURPOSE:
C WLOG3D: Compute stress and tangent in principal stretch form
C
C METHOD:
C Stored Energy Function in Principal Stretches:
C W = w(lam_1) + w(lam_2) + w(lam_3)
C taup_i = [w(lam_i)]' * lambda_i
C gamm_ij = lambda_j * [taup_i],j
C
C Logarithmic strain energy function
C w = mu * [log(lambda)]**2
C
C ARGUMENTS INPUT:
C mu - Shear modulus
C bpr(3) - Principal stretch (squared)
C
C ARGUMENTS OUTPUT:
C epsd(3) - Principal deviatoric strains
C taup(3) - Principal deviatoric Kirchhoff stresses
C aap(6,6) - Deviatoric tangent moduli in principal basis
C ierr - Error flag
C = 1 :
C
C COMMON BLOCKS:
C const
C
C PERIPHERAL UNITS:
C None
C
C INTERNAL VARIABLES:
C
C
C
C PRINT SWITCH:
C ipswb3 = 0 Gives no print
C ipswb3 = 2 Gives enter and leave
C ipswb3 = 3 Gives in addition parameters on input
C ipswb3 = 5 Gives in addition parameters on output
C
C LIMITS:
C
C CODING:
C STANDART ADHERED TO: Computas report no 78-949
C AUTHOR: Kjell Magne Mathisen
C NTNU, Department of structural engineering
C November 14, 2010
C
C INSPECTED:
C
C REVISIONS:
C xx-OCT-2010 K.M.Mathisen NTNU
C Changed XG1-XG8 to XG and UG1-UG8 to UG.
C
C ---------------------------------------------------------------------
C
implicit none
C
integer i, j, ipswb3, iwr, ierr
C
real*8 jthird, mu, twomu, fourmu, w, tol,
& bpr(3), lamt(3), epsd(3), taup(3), aap(6,6)
C
include 'include/feninc/const.h'
C
data tol / 1.0d-08 /
C
C Entry section
C
ierr = 0
C
if (ipswb3 .gt. 0) then
write(iwr,9010) 'ENTERING SUBROUTINE WLOG3D'
if (ipswb3 .gt. 2) then
write(iwr,9010) 'WITH INPUT ARGUMENTS'
write(iwr,9030) 'mu =', mu
call rprin0(bpr, 1, 3, 'bpr ', iwr)
endif
endif
C
C Compute Je**(-1/3) and volume preserving stretches
C
jthird = (bpr(1)*bpr(2)*bpr(3))**(-one6)
do i = 1,3
lamt(i) = jthird*sqrt(bpr(i))
end do ! i
C
C Compute predictor strain energy function derivatives
C
twomu = two * mu
w = zero
do i = 1,3
epsd(i) = log(lamt(i))
taup(i) = twomu * epsd(i)
w = w + mu * epsd(i)**2
end do ! i
C
C Load material part of moduli in principal stretches
C
aap = zero
C
twomu = twomu * one3
fourmu = twomu + twomu
C
do i = 1,3
aap(i,i) = fourmu - two * taup(i)
end do ! i
C
aap(1,2) = - twomu
aap(2,1) = - twomu
C
aap(2,3) = - twomu
aap(3,2) = - twomu
C
aap(3,1) = - twomu
aap(1,3) = - twomu
C
C Compute deviatoric tangent moduli in principal basis
C
do i = 1,3
j = mod(i,3) + 1
C
if (abs(bpr(i) - bpr(j)) .gt. tol) then
aap(i+3,i+3) = ( bpr(i)*taup(j) - bpr(j)*taup(i) )
& / ( bpr(j) - bpr(i) )
else
aap(i+3,i+3) = (aap(i,i) - aap(j,i))*one2
endif
C
end do ! i
go to 8000
C
C Error section
C
7000 continue
ierr =-1
go to 7900
7100 continue
ierr =-2
go to 7900
C ----------
7900 continue
write(iwr,9020) '*** ERROR IN SUBROUTINE WLOG3D *** IERR = ', ierr
write(iwr,9010) 'WITH INPUT ARGUMENTS'
write(iwr,9030) 'mu =', mu
call rprin0(bpr, 1, 3, 'bpr ', iwr)
C
write(iwr,9010) 'WITH OUTPUT ARGUMENTS'
write(iwr,9030) 'w =', w
call rprin0(epsd, 1, 3, 'epsd ', iwr)
call rprin0(taup, 1, 3, 'taup ', iwr)
call rprin0(aap , 6, 6, 'aap ', iwr)
go to 8010
C
C Closing section
C
8000 continue
C
if (ipswb3 .gt. 0) then
write(iwr,9010) 'LEAVING SUBROUTINE WLOG3D'
if (ipswb3 .gt. 3) then
write(iwr,9010) 'WITH OUTPUT ARGUMENTS'
write(iwr,9030) 'w =', w
call rprin0(epsd, 1, 3, 'epsd ', iwr)
call rprin0(taup, 1, 3, 'taup ', iwr)
call rprin0(aap , 6, 6, 'aap ', iwr)
endif
endif
C
8010 continue
C
return
C
9010 format(3X,A)
9020 format(3X,A,I6)
9030 format(3X,A,E15.6)
C
end

View File

@@ -0,0 +1,187 @@
// $Id: NeoHookeElasticity.C,v 1.2 2011-02-08 09:06:02 kmo Exp $
//==============================================================================
//!
//! \file NeoHookeElasticity.C
//!
//! \date Jul 1 2010
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Integrand implementations for Neo-Hooke elasticity problems.
//!
//==============================================================================
#include "NeoHookeElasticity.h"
void NeoHookeElasticity::print (std::ostream& os) const
{
this->NonlinearElasticity::print(os);
std::cout <<"NeoHookeElasticity: Neo-Hooke material formulation"<< std::endl;
}
void NeoHookeElasticity::setMaterial (double Emod, double Poiss, double Density)
{
this->Elasticity::setMaterial(Emod,Poiss,Density);
// Define the Lame elasticity parameters
lambda = Emod*Poiss / ((1.0+Poiss)*(1.0-Poiss-Poiss));
mu = lambda * (0.5/Poiss-1.0);
kappa = lambda + (mu+mu)/3.0;
}
bool NeoHookeElasticity::formStressTensor (const Matrix& dNdX, const Vec3&,
SymmTensor& S) const
{
if (!eV || eV->empty())
{
// Initial state (zero stresses)
S.zero();
return true;
}
// Evaluate the kinematic quantities, F, E, Ci, and J at this point
if (!this->kinematics(dNdX,E))
return false;
// Set up the 2nd Piola-Kirchhoff stress tensor, S
this->formStressTensor(S);
return true;
}
bool NeoHookeElasticity::formTangent (Matrix& Ctan, SymmTensor& S,
const Vec3& X) const
{
if (!this->formCmatrix(Ctan,X))
return false;
else if (eV && !eV->empty())
this->formStressTensor(S);
return true;
}
bool NeoHookeElasticity::kinematics (const Matrix& dNdX, SymmTensor& _E) const
{
// Form the deformation gradient, F
if (!this->NonlinearElasticity::kinematics(dNdX,_E))
return false;
else
J = F.det();
// Form the right Cauchy-Green tensor, C = F^T*F
static Tensor dUdX(nsd);
C.rightCauchyGreen(dUdX=F);
// Invert the right Cauchy-Green tensor, Ci = C^-1
Ci = C;
return Ci.inverse() > 0.0;
}
void NeoHookeElasticity::formStressTensor (SymmTensor& S) const
{
// Set up the 2nd Piola-Kirchhoff stress tensor.
// See equation (3.119) in P. Wriggers book:
// S = mu*(1-Ci) + 0.5*lambda*(J^2-1)*Ci
S = 0.5*lambda*(J*J-1.0)*Ci;
S += -mu*(Ci-1.0);
}
bool NeoHookeElasticity::formCmatrix (Matrix& Ctan, const Vec3&, bool) const
{
// Set up the incremental consitutive tensor
size_t nst = nsd*(nsd+1)/2;
Ctan.resize(nst,nst,true);
SymmTensor4 D(Ctan,nsd); // fourth-order material tensor
double a = lambda*J*J;
double b = mu - 0.5*(a-lambda);
unsigned short int i, j, k, l;
for (i = 1; i <= nsd; i++)
for (j = i; j <= nsd; j++)
for (k = 1; k <= nsd; k++)
for (l = k; l <= nsd; l++)
// See equation (3.268) in P. Wriggers book
D(i,j,k,l) = a*Ci(i,j)*Ci(k,l) + b*(Ci(i,k)*Ci(j,l) + Ci(i,l)*C(j,k));
return true;
}
void NeoHookeElasticityIV::print (std::ostream& os) const
{
this->NonlinearElasticity::print(os);
std::cout <<"NeoHookeElasticity: Neo-Hooke material formulation"
<<" with isochoric/volumetric decomposition"<< std::endl;
}
bool NeoHookeElasticityIV::formTangent (Matrix& Ctan, SymmTensor& S,
const Vec3& X) const
{
// Set up the stress tensor first because Ctan also depends on Siso
if (eV && !eV->empty())
this->formStressTensor(S);
else
{
// Initial state (zero stresses)
S.zero();
Siso.zero();
}
return this->formCmatrix(Ctan,X);
}
void NeoHookeElasticityIV::formStressTensor (SymmTensor& S) const
{
// Set up the 2nd Piola-Kirchhoff stress tensor.
// See equation (3.127) or (3.274) in P. Wriggers book:
// S = mu*J^(-2/3)*(1-tr(C)/3)*Ci
S = -(C.trace()/3.0)*Ci;
S += 1.0;
S *= mu*pow(J,-2.0/3.0);
Siso = S; // Isochoric part
S += 0.5*kappa*(J*J-1.0)*Ci; // Volumetric part
}
bool NeoHookeElasticityIV::formCmatrix (Matrix& Ctan, const Vec3&, bool) const
{
// Set up the incremental consitutive tensor
size_t nst = nsd*(nsd+1)/2;
Ctan.resize(nst,nst,true);
SymmTensor4 D(Ctan,nsd); // fourth-order material tensor
double t = mu*pow(J,-2.0/3.0)*C.trace();
double a = kappa*J*J + (t+t)/9.0;
double b = 0.5*kappa*(1.0-J*J) + t/3.0;
double c = (t+t)/3.0; // Yuri
//double c = 2.0/3.0; // Wriggers
unsigned short int i, j, k, l;
for (i = 1; i <= nsd; i++)
for (j = i; j <= nsd; j++)
for (k = 1; k <= nsd; k++)
for (l = k; l <= nsd; l++)
{
// See equation (3.275) in P. Wriggers book
//D(i,j,k,l) = a*Ci(i,j)*Ci(k,l) + b*(Ci(i,k)*Ci(j,l) + Ci(i,l)*C(j,k))
// - c*(Ci(k,l)*Siso(i,j) + Ci(i,j)*Siso(k,l));
// See equation (81) In Yuri's Comput. Mech. paper
D(i,j,k,l) = a*Ci(i,j)*Ci(k,l) + b*(Ci(i,k)*Ci(j,l) + Ci(i,l)*C(j,k));
if (i == j) D(i,j,k,l) += c*Ci(k,l);
if (k == l) D(i,j,k,l) += c*Ci(i,j);
}
return true;
}

View File

@@ -0,0 +1,128 @@
// $Id: NeoHookeElasticity.h,v 1.1 2010-12-18 20:03:53 kmo Exp $
//==============================================================================
//!
//! \file NeoHookeElasticity.h
//!
//! \date Jul 1 2010
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Integrand implementations for Neo-Hooke elasticity problems.
//!
//==============================================================================
#ifndef _NEO_HOOKE_ELASTICITY_H
#define _NEO_HOOKE_ELASTICITY_H
#include "NonlinearElasticity.h"
/*!
\brief Class representing the integrand of nonlinear hyperelastic problems.
\details This class implements the Neo-Hooke hyperelastic material behavior.
*/
class NeoHookeElasticity : public NonlinearElasticity
{
public:
//! \brief The default constructor invokes the parent class constructor only.
//! \param[in] n Number of spatial dimensions
NeoHookeElasticity(unsigned short int n = 3)
: NonlinearElasticity(n), C(n), Ci(n) { fullCmat = true; }
//! \brief Empty destructor.
virtual ~NeoHookeElasticity() {}
//! \brief Prints out problem definition to the given output stream.
virtual void print(std::ostream& os) const;
//! \brief Defines material properties for current volume patch.
//! \param[in] Emod Young's modulus
//! \param[in] Poiss Poisson's ratio
//! \param[in] Density Mass density
virtual void setMaterial(double Emod, double Poiss, double Density);
protected:
//! \brief Calculates some kinematic quantities at current point.
//! \param[in] dNdX Basis function gradients at current point
//! \param[out] _E Green-Lagrange strain tensor at current point
//!
//! \details The deformation gradient \b F and the Green-Lagrange tensor \b E
//! are established. In addition, the right Cauchy-Green tensor \b C and its
//! inverse \b Ci are formed.
//! All quantities are stored internally in mutable data members.
virtual bool kinematics(const Matrix& dNdX, SymmTensor& _E) const;
//! \brief Forms tangential tensorial quantities needed by the evalInt method.
//! \param[in] X Cartesian coordinates of current integration point
//! \param[out] Ctan Tangential constitutive tensor at current point
//! \param[out] S 2nd Piola-Kirchhoff stress tensor at current point
virtual bool formTangent(Matrix& Ctan, SymmTensor& S, const Vec3& X) const;
//! \brief Forms the 2nd Piola-Kirchhoff stress tensor.
//! \param[in] dNdX Basis function gradients at current integration point
//! \param[in] X Cartesian coordinates of current integration point
//! \param[out] S 2nd Piola-Kirchhoff stress tensor at current point
virtual bool formStressTensor(const Matrix& dNdX, const Vec3& X,
SymmTensor& S) const;
//! \brief Forms the 2nd Piola-Kirchhoff stress tensor.
virtual void formStressTensor(SymmTensor& S) const;
public:
//! \brief Sets up the tangential constitutive matrix at current point.
//! \param[out] Ctan \f$6\times6\f$-matrix (in 3D) or \f$3\times3\f$-matrix
//! (in 2D), representing the tangential constitutive tensor
//! \param[in] X Cartesian coordinates of current point
virtual bool formCmatrix(Matrix& Ctan, const Vec3& X, bool = false) const;
protected:
double lambda; //!< First Lame parameter
double mu; //!< Second Lame parameter (shear modulus)
double kappa; //!< Bulk modulus
mutable double J; //!< Determinant of the deformation gradient
mutable SymmTensor C; //!< Right Cauchy-Green deformation tensor
mutable SymmTensor Ci; //!< Inverse of the right Cauchy-Green tensor
};
/*!
\brief Neo-Hooke hyperelastic material with isochoric and volumetric split.
*/
class NeoHookeElasticityIV : public NeoHookeElasticity
{
public:
//! \brief The default constructor invokes the parent class constructor only.
//! \param[in] n Number of spatial dimensions
NeoHookeElasticityIV(unsigned short int n = 3)
: NeoHookeElasticity(n), Siso(n) {}
//! \brief Empty destructor.
virtual ~NeoHookeElasticityIV() {}
//! \brief Prints out problem definition to the given output stream.
virtual void print(std::ostream& os) const;
protected:
//! \brief Forms tangential tensorial quantities needed by the evalInt method.
//! \param[in] X Cartesian coordinates of current integration point
//! \param[out] Ctan Tangential constitutive tensor at current point
//! \param[out] S 2nd Piola-Kirchhoff stress tensor at current point
virtual bool formTangent(Matrix& Ctan, SymmTensor& S, const Vec3& X) const;
//! \brief Forms the 2nd Piola-Kirchhoff stress tensor.
virtual void formStressTensor(SymmTensor& S) const;
public:
//! \brief Sets up the tangential constitutive matrix at current point.
//! \param[out] Ctan \f$6\times6\f$-matrix (in 3D) or \f$3\times3\f$-matrix
//! (in 2D), representing the tangential constitutive tensor
//! \param[in] X Cartesian coordinates of current point
virtual bool formCmatrix(Matrix& Ctan, const Vec3& X, bool = false) const;
protected:
mutable SymmTensor Siso; //!< Isochoric part of the 2nd PK stress tensor
};
#endif

View File

@@ -0,0 +1,250 @@
// $Id: NonlinearElasticity.C,v 1.3 2011-02-08 09:06:02 kmo Exp $
//==============================================================================
//!
//! \file NonlinearElasticity.C
//!
//! \date May 25 2010
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Integrand implementations for nonlinear elasticity problems.
//!
//==============================================================================
#include "NonlinearElasticity.h"
#include "Utilities.h"
#include "Profiler.h"
#ifndef NO_FORTRAN
#if defined(_WIN32)
#define stiff_tl2d_ STIFF_TL2D
#define stiff_tl3d_ STIFF_TL3D
#define stiff_tl3d_isoel_ STIFF_TL3D_ISOEL
#define stiff_tl3d_isoel_ STIFF_TL3D_ISOEL
#elif defined(_AIX)
#define stiff_tl2d_ stiff_tl2d
#define stiff_tl3d_ stiff_tl3d
#define stiff_tl2d_isoel_ stiff_tl2d_isoel
#define stiff_tl3d_isoel_ stiff_tl3d_isoel
#endif
extern "C" {
//! \brief Calculates material stiffness contributions for 2D problems.
void stiff_tl2d_(const int& nenod, const double& detJW,
const double* dNdX, const double* F, const double* Cmat,
double* EM);
//! \brief Calculates material stiffness contributions for 3D problems.
void stiff_tl3d_(const int& nenod, const double& detJW,
const double* dNdX, const double* F, const double* Cmat,
double* EM);
//! \brief Optimized for isotropic linear elastic materials in 2D.
void stiff_tl2d_isoel_(const int& nenod, const double& detJW,
const double* dNdX, const double* F,
const double& C1, const double& C2,const double& C3,
double* EM);
//! \brief Optimized for isotropic linear elastic materials in 3D.
void stiff_tl3d_isoel_(const int& nenod, const double& detJW,
const double* dNdX, const double* F,
const double& C1, const double& C2,const double& C3,
double* EM);
}
#endif
NonlinearElasticity::NonlinearElasticity (unsigned short int n, bool ps)
: NonlinearElasticityTL(n,ps), E(n)
{
fullCmat = false;
}
void NonlinearElasticity::print (std::ostream& os) const
{
this->Elasticity::print(os);
std::cout <<"NonlinearElasticity: Total Lagrangian formulation "
<<" (tensorial form)"<< std::endl;
}
void NonlinearElasticity::setMode (SIM::SolutionMode mode)
{
this->NonlinearElasticityTL::setMode(mode);
formB = false; // We don't need the B-matrix in the tensor formulation
}
bool NonlinearElasticity::evalInt (LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const
{
PROFILE3("NonlinearEl::evalInt");
// Evaluate the kinematic quantities, F and E, at this point
if (!this->kinematics(dNdX,E))
return false;
// Evaluate current tangent at this point, that is
// the incremental constitutive matrix, Cmat, and
// the 2nd Piola-Kirchhoff stress tensor, S
static SymmTensor S(nsd);
if (!this->formTangent(Cmat,S,X))
return false;
bool haveStrains = !E.isZero(1.0e-16);
size_t a, b;
unsigned short int i, j, k;
if (eKm)
{
// Integrate the material stiffness matrix
#ifndef NO_FORTRAN
// Using Fortran routines optimized for symmetric constitutive tensors
PROFILE4("stiff_TL_");
if (nsd == 3)
if (fullCmat)
stiff_tl3d_(N.size(),detJW,dNdX.ptr(),F.ptr(),
Cmat.ptr(),eKm->ptr());
else
stiff_tl3d_isoel_(N.size(),detJW,dNdX.ptr(),F.ptr(),
Cmat(1,1),Cmat(1,2),Cmat(4,4),eKm->ptr());
else if (nsd == 2)
if (fullCmat)
stiff_tl2d_(N.size(),detJW,dNdX.ptr(),F.ptr(),
Cmat.ptr(),eKm->ptr());
else
stiff_tl2d_isoel_(N.size(),detJW,dNdX.ptr(),F.ptr(),
Cmat(1,1),Cmat(1,2),Cmat(3,3),eKm->ptr());
else if (nsd == 1)
for (a = 1; a <= N.size(); a++)
for (b = 1; b <= N.size(); b++)
(*eKm)(a,b) += dNdX(a,1)*F(1,1)*Cmat(1,1)*F(1,1)*dNdX(b,1);
#else
// This is too costly, but is basically what is done in the fortran routines
PROFILE4("dNdX^t*F^t*C*F*dNdX");
unsigned short int l, m, n;
SymmTensor4 D(Cmat,nsd); // fourth-order material tensor
Matrix& EM = *eKm;
for (a = 1; a <= N.size(); a++)
for (b = 1; b <= N.size(); b++)
for (m = 1; m <= nsd; m++)
for (n = 1; n <= nsd; n++)
{
double km = 0.0;
for (i = 1; i <= nsd; i++)
for (j = 1; j <= nsd; j++)
for (k = 1; k <= nsd; k++)
for (l = 1; l <= nsd; l++)
km += dNdX(a,i)*F(m,j)*D(i,j,k,l)*F(n,k)*dNdX(b,l);
EM(nsd*(a-1)+m,nsd*(b-1)+n) += km*detJW;
}
#endif
}
if (eKg && haveStrains)
// Integrate the geometric stiffness matrix
this->formKG(*eKg,dNdX,S,detJW);
if (eM)
// Integrate the mass matrix
this->formMassMatrix(*eM,N,X,detJW);
if (iS && haveStrains)
{
// Integrate the internal forces
Vector& ES = *iS;
for (a = 1; a <= N.size(); a++)
for (k = 1; k <= nsd; k++)
{
double f = 0.0;
for (i = 1; i <= nsd; i++)
for (j = 1; j <= nsd; j++)
f -= dNdX(a,i)*F(k,j)*S(i,j);
ES(nsd*(a-1)+k) += f*detJW;
}
}
if (eS)
// Integrate the load vector due to gravitation and other body forces
this->formBodyForce(*eS,N,X,detJW);
return this->getIntegralResult(elmInt);
}
bool NonlinearElasticity::evalSol (Vector& s, const Vector&,
const Matrix& dNdX, const Vec3& X,
const std::vector<int>& MNPC) const
{
PROFILE3("NonlinearEl::evalSol");
int ierr = 0;
if (eV && !primsol.front().empty())
if (ierr = utl::gather(MNPC,nsd,primsol.front(),*eV))
{
std::cerr <<" *** NonlinearElasticity::evalSol: Detected "
<< ierr <<" node numbers out of range."<< std::endl;
return false;
}
// Evaluate the stress state at this point
static SymmTensor Sigma(nsd);
if (!this->formStressTensor(dNdX,X,Sigma))
return false;
// Congruence transformation to local coordinate system at current point
if (locSys) Sigma.transform(locSys->getTmat(X));
s = Sigma;
s.push_back(Sigma.vonMises());
return true;
}
bool NonlinearElasticity::evalSol (Vector& s,
const Matrix& dNdX, const Vec3& X) const
{
PROFILE3("NonlinearEl::evalSol");
static SymmTensor Sigma(nsd);
if (!this->formStressTensor(dNdX,X,Sigma))
return false;
s = Sigma;
return true;
}
bool NonlinearElasticity::formStressTensor (const Matrix& dNdX, const Vec3& X,
SymmTensor& S) const
{
if (!eV || eV->empty())
{
// Initial state (zero stresses)
S.zero();
return true;
}
// Evaluate the kinematic quantities, F and E, at this point
if (!this->kinematics(dNdX,E))
return false;
// Evaluate the constitutive matrix, C, at this point
if (!this->formCmatrix(Cmat,X))
return false;
return Cmat.multiply(E,S); // S = C*E
}
bool NonlinearElasticity::formTangent (Matrix& Ctan, SymmTensor& S,
const Vec3& X) const
{
if (!this->formCmatrix(Ctan,X))
return false;
else if (eV && !eV->empty())
return Ctan.multiply(E,S); // S = C*E
else
return true; // Initial state (no stresses)
}

View File

@@ -0,0 +1,96 @@
// $Id: NonlinearElasticity.h,v 1.2 2011-02-08 09:06:02 kmo Exp $
//==============================================================================
//!
//! \file NonlinearElasticity.h
//!
//! \date May 25 2010
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Integrand implementations for nonlinear elasticity problems.
//!
//==============================================================================
#ifndef _NONLINEAR_ELASTICITY_H
#define _NONLINEAR_ELASTICITY_H
#include "NonlinearElasticityTL.h"
#include "Tensor.h"
/*!
\brief Class representing the integrand of the nonlinear elasticity problem.
\details This class implements a Total Lagrangian formulation, tensorial form.
It reimplements most of the Elasticity methods, except for the methods
\a kinematics and \a evalBou, which are inherited from NonlinearElasticityTL.
\note This class is obsolete, as it gives exactly the same results as the
NonlinearElasticityTL class, but is less efficient. It is retained mostly for
educational purposes and historical reasons, and such that the performance of
the two ways of implementing the total Lagrangian formulation can be compared
on a varity of problems.
*/
class NonlinearElasticity : public NonlinearElasticityTL
{
public:
//! \brief The default constructor invokes the parent class constructor only.
//! \param[in] n Number of spatial dimensions
//! \param[in] ps If \e true, assume plane stress in 2D
NonlinearElasticity(unsigned short int n = 3, bool ps = true);
//! \brief Empty destructor.
virtual ~NonlinearElasticity() {}
//! \brief Prints out problem definition to the given output stream.
virtual void print(std::ostream& os) const;
//! \brief Defines the solution mode before the element assembly is started.
//! \param[in] mode The solution mode to use
virtual void setMode(SIM::SolutionMode mode);
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] X Cartesian coordinates of current integration point
virtual bool evalInt(LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const;
//! \brief Evaluates the secondary solution at a result point.
//! \param[out] s Array of solution field values at current point
//! \param[in] N Basis function values at current point
//! \param[in] dNdX Basis function gradients at current point
//! \param[in] X Cartesian coordinates of current point
//! \param[in] MNPC Nodal point correspondance for the basis function values
virtual bool evalSol(Vector& s,
const Vector& N, const Matrix& dNdX,
const Vec3& X, const std::vector<int>& MNPC) const;
//! \brief Evaluates the finite element (FE) solution at an integration point.
//! \param[out] s The FE stress values at current point
//! \param[in] dNdX Basis function gradients at current point
//! \param[in] X Cartesian coordinates of current point
virtual bool evalSol(Vector& s, const Matrix& dNdX, const Vec3& X) const;
protected:
//! \brief Forms tangential tensorial quantities needed by the evalInt method.
//! \param[in] X Cartesian coordinates of current integration point
//! \param[out] Ctan Tangential constitutive tensor at current point
//! \param[out] S 2nd Piola-Kirchhoff stress tensor at current point
virtual bool formTangent(Matrix& Ctan, SymmTensor& S, const Vec3& X) const;
//! \brief Forms the 2nd Piola-Kirchhoff stress tensor.
//! \param[in] dNdX Basis function gradients at current integration point
//! \param[in] X Cartesian coordinates of current integration point
//! \param[out] S 2nd Piola-Kirchhoff stress tensor at current point
virtual bool formStressTensor(const Matrix& dNdX, const Vec3& X,
SymmTensor& S) const;
protected:
bool fullCmat; //!< If \e true, assume a full (but symmetric) C-matrix
mutable SymmTensor E; //!< Green-Lagrange strain tensor
};
#endif

View File

@@ -0,0 +1,221 @@
// $Id: NonlinearElasticityTL.C,v 1.4 2011-02-08 09:06:02 kmo Exp $
//==============================================================================
//!
//! \file NonlinearElasticityTL.C
//!
//! \date May 25 2010
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Integrand implementations for nonlinear elasticity problems.
//!
//==============================================================================
#include "NonlinearElasticity.h"
#include "ElmMats.h"
#include "Tensor.h"
#include "Vec3Oper.h"
NonlinearElasticityTL::NonlinearElasticityTL (unsigned short int n, bool ps)
: Elasticity(n,ps)
{
// Only the current solution is needed
primsol.resize(1);
}
void NonlinearElasticityTL::print (std::ostream& os) const
{
this->Elasticity::print(os);
std::cout <<"NonlinearElasticityTL: Total Lagrangian formulation"<< std::endl;
}
void NonlinearElasticityTL::setMode (SIM::SolutionMode mode)
{
if (!myMats) return;
formB = true;
myMats->rhsOnly = false;
eM = eKm = eKg = 0;
iS = eS = eV = 0;
switch (mode)
{
case SIM::STATIC:
myMats->resize(1,2);
eKm = &myMats->A[0];
eKg = &myMats->A[0];
break;
case SIM::DYNAMIC:
myMats->resize(2,2);
eKm = &myMats->A[0];
eKg = &myMats->A[0];
eM = &myMats->A[1];
break;
case SIM::RHS_ONLY:
if (myMats->A.empty())
myMats->resize(1,2);
else
myMats->b.resize(2);
eKm = &myMats->A[0];
eKg = &myMats->A[0];
myMats->rhsOnly = true;
break;
default:
formB = false;
this->Elasticity::setMode(mode);
return;
}
// We always need the force+displacement vectors in nonlinear simulations
iS = &myMats->b[0];
eS = &myMats->b[0];
eV = &myMats->b[1];
tracVal.clear();
}
bool NonlinearElasticityTL::evalBou (LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X, const Vec3& normal) const
{
if (!tracFld)
{
std::cerr <<" *** NonlinearElasticityTL::evalBou: No tractions."
<< std::endl;
return false;
}
else if (!eS)
{
std::cerr <<" *** NonlinearElasticityTL::evalBou: No load vector."
<< std::endl;
return false;
}
// Evaluate the surface traction
Vec3 T = (*tracFld)(X,normal);
// Store the traction value for vizualization
if (!T.isZero()) tracVal[X] = T;
// Check for with-rotated pressure load
unsigned short int i, j;
if (tracFld->isNormalPressure())
{
// Compute the deformation gradient, F
SymmTensor dummy(0);
if (!this->kinematics(dNdX,dummy)) return false;
// Compute its inverse and determinant, J
double J = F.inverse();
if (J <= 0.0) return false;
// Pull-back the normal traction to the initial configuration.
// See equation (3.4.5) on page 102 in T. Belytschko et. al (2000):
// p*n*dS ==> J*p*n0*F^-1*dS0
Vec3 t = J*T; T = 0.0;
for (i = 1; i <= nsd; i++)
for (j = 1; j <= nsd; j++)
T[i-1] += t[j-1]*F(j,i);
}
// Integrate the force vector
Vector& ES = *eS;
for (size_t a = 1; a <= N.size(); a++)
for (i = 1; i <= nsd; i++)
ES(nsd*(a-1)+i) += T[i-1]*N(a)*detJW;
return this->getIntegralResult(elmInt);
}
bool NonlinearElasticityTL::kinematics (const Matrix& dNdX, SymmTensor& E) const
{
if (!eV || eV->empty())
{
// Initial state, unit deformation gradient and linear B-matrix
F.diag(1.0,nsd);
return this->Elasticity::kinematics(dNdX,E);
}
const size_t nenod = dNdX.rows();
const size_t nstrc = nsd*(nsd+1)/2;
if (eV->size() != nenod*nsd || dNdX.cols() < nsd)
{
std::cerr <<" *** NonlinearElasticityTL::kinematics: Invalid dimension,"
<<" dNdX("<< nenod <<","<< dNdX.cols() <<")"<< std::endl;
return false;
}
// Compute the deformation gradient, [F] = [I] + [dudX] = [I] + [dNdX]*[u],
// and the Green-Lagrange strains, E_ij = 0.5(F_ij+F_ji+F_ki*F_kj).
// Notice that the matrix multiplication method used here treats the
// element displacement vector, *eV, as a matrix whose number of columns
// equals the number of rows in the matrix dNdX.
if (!F.multiplyMat(*eV,dNdX)) // F = Grad{u} = eV*dNdX
return false;
unsigned short int i, j, k;
// Now form the Green-Lagrange strain tensor.
// Note that for the shear terms (i/=j) we actually compute 2*E_ij
// to be consistent with the engineering strain style constitutive matrix.
for (i = 1; i <= E.dim(); i++)
for (j = 1; j <= i; j++)
{
double Eij = F(i,j) + F(j,i);
for (k = 1; k <= nsd; k++)
Eij += F(k,i)*F(k,j);
E(i,j) = i == j ? 0.5*Eij : Eij;
}
// Add the unit tensor to F to form the deformation gradient
for (i = 1; i <= nsd; i++)
F(i,i) += 1.0;
#ifdef INT_DEBUG
std::cout <<"NonlinearElasticityTL::F ="<< F;
#endif
if (!formB || E.dim() < nsd) return true;
// Form the nonlinear B-matrix
Bmat.resize(nstrc*nsd,nenod,true);
#define INDEX(i,j) i+nstrc*(j-1)
// Normal strain part
size_t a;
for (a = 1; a <= nenod; a++)
for (i = 1; i <= nsd; i++)
for (j = 1; j <= nsd; j++)
Bmat(INDEX(j,i),a) = F(i,j)*dNdX(a,j);
// Shear strain part
if (nsd == 3)
for (a = 1; a <= nenod; a++)
for (i = 1; i <= nsd; i++)
{
Bmat(INDEX(4,i),a) = F(i,1)*dNdX(a,2) + F(i,2)*dNdX(a,1);
Bmat(INDEX(5,i),a) = F(i,2)*dNdX(a,3) + F(i,3)*dNdX(a,2);
Bmat(INDEX(6,i),a) = F(i,3)*dNdX(a,1) + F(i,1)*dNdX(a,3);
}
else if (nsd == 2)
for (a = 1; a <= nenod; a++)
for (i = 1; i <= nsd; i++)
Bmat(INDEX(3,i),a) = F(i,1)*dNdX(a,2) + F(i,2)*dNdX(a,1);
Bmat.resize(nstrc,nsd*nenod);
#ifdef INT_DEBUG
std::cout <<"NonlinearElasticityTL::B ="<< Bmat;
#endif
return true;
}

View File

@@ -0,0 +1,76 @@
// $Id: NonlinearElasticityTL.h,v 1.2 2011-02-08 09:06:02 kmo Exp $
//==============================================================================
//!
//! \file NonlinearElasticityTL.h
//!
//! \date May 25 2010
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Integrand implementations for nonlinear elasticity problems.
//!
//==============================================================================
#ifndef _NONLINEAR_ELASTICITY_TL_H
#define _NONLINEAR_ELASTICITY_TL_H
#include "Elasticity.h"
/*!
\brief Class representing the integrand of the nonlinear elasticity problem.
\details This class implements a Total Lagrangian formulation in matrix form.
It inherits most of the Elasticity methods, but reimplements the \a kinematics
method, for calculating the nonlinear variant of the strain-displacement
matrix, \b B, and the associated Green-Lagrange strain tensor, \b E.
The \a evalBou method is also reimplemented to account for with-rotated loads.
*/
class NonlinearElasticityTL : public Elasticity
{
public:
//! \brief The default constructor invokes the parent class constructor only.
//! \param[in] n Number of spatial dimensions
//! \param[in] ps If \e true, assume plane stress in 2D
NonlinearElasticityTL(unsigned short int n = 3, bool ps = true);
//! \brief Empty destructor.
virtual ~NonlinearElasticityTL() {}
//! \brief Prints out problem definition to the given output stream.
virtual void print(std::ostream& os) const;
//! \brief Defines the solution mode before the element assembly is started.
//! \param[in] mode The solution mode to use
virtual void setMode(SIM::SolutionMode mode);
//! \brief Evaluates the integrand at a boundary point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] X Cartesian coordinates of current integration point
//! \param[in] normal Boundary normal vector at current integration point
//!
//! \details This method is reimplemented in this class to account for
//! possibly with-rotated traction fields in the Total-Lagrangian setting.
virtual bool evalBou(LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X, const Vec3& normal) const;
protected:
//! \brief Calculates some kinematic quantities at current point.
//! \param[in] dNdX Basis function gradients at current point
//! \param[out] E Green-Lagrange strain tensor at current point
//!
//! \details The deformation gradient \b F and the nonlinear
//! strain-displacement matrix \b B are established. The latter matrix
//! is stored in the mutable class member \a Bmat of the parent class.
//! The B-matrix is formed only when the variable \a formB is true.
virtual bool kinematics(const Matrix& dNdX, SymmTensor& E) const;
protected:
bool formB; //!< Flag determining whether we need to form the B-matrix
mutable Matrix F; //!< Deformation gradient
};
#endif

View File

@@ -0,0 +1,468 @@
// $Id: NonlinearElasticityUL.C,v 1.5 2011-02-08 09:06:02 kmo Exp $
//==============================================================================
//!
//! \file NonlinearElasticityUL.C
//!
//! \date Sep 21 2010
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Integrand implementations for nonlinear elasticity problems.
//!
//==============================================================================
#include "NonlinearElasticityUL.h"
#include "ElmMats.h"
#include "ElmNorm.h"
#include "Tensor.h"
#include "Vec3Oper.h"
#ifdef USE_FTNMAT
extern "C" {
//! \brief Interface to 2D nonlinear material routines (FORTRAN-77 code).
void cons2d_(const int& mTYP, const int& mVER, const double& detF,
const double* F, const double* pmat, double& Engy,
const double* Sig, double* Cst,
const int& ipsw, const int& iwr, int& ierr);
//! \brief Interface to 3D nonlinear material routines (FORTRAN-77 code).
void cons3d_(const int& mTYP, const int& mVER, const double& detF,
const double* F, const double* pmat, double& Engy,
const double* Sig, double* Cst,
const int& ipsw, const int& iwr, int& ierr);
}
#ifndef INT_DEBUG
#define INT_DEBUG 0
#endif
#endif
NonlinearElasticityUL::NonlinearElasticityUL (unsigned short int n,
bool ps, int mver, char lop)
: Elasticity(n,ps)
{
mVER = -1;
#ifdef USE_FTNMAT
if (n == 3 || (n == 2 && !ps)) mVER = mver;
#endif
loadOp = lop;
// Only the current solution is needed
primsol.resize(1);
}
void NonlinearElasticityUL::print (std::ostream& os) const
{
this->Elasticity::print(os);
std::cout <<"NonlinearElasticityUL: Updated Lagrangian formulation";
if (mVER >= 0) std::cout <<", mVER="<< mVER;
if (mVER >= 1) std::cout <<" (Neo-Hooke)";
std::cout << std::endl;
}
void NonlinearElasticityUL::setMaterial (double Young, double Poiss,
double Density)
{
mTYP = 10;
pmat[0] = Young;
pmat[1] = Poiss;
if (Poiss > 0.5)
{
mTYP *= -1; // Assume Lame' parameters (kappa and mu are specified)
double kappa = Young;
double mu = Poiss;
// Calculate Young's modulus and Poisson's ratio
Young = 9.0*kappa*mu/(3.0*kappa + mu);
Poiss = (1.5*kappa - mu)/(3.0*kappa + mu);
}
this->Elasticity::setMaterial(Young,Poiss,Density);
}
void NonlinearElasticityUL::setMode (SIM::SolutionMode mode)
{
if (!myMats) return;
myMats->rhsOnly = false;
eM = eKm = eKg = 0;
iS = eS = eV = 0;
switch (mode)
{
case SIM::STATIC:
myMats->resize(1,2);
eKm = &myMats->A[0];
eKg = &myMats->A[0];
break;
case SIM::DYNAMIC:
myMats->resize(2,2);
eKm = &myMats->A[0];
eKg = &myMats->A[0];
eM = &myMats->A[1];
break;
case SIM::RHS_ONLY:
if (myMats->A.empty())
myMats->resize(1,2);
else
myMats->b.resize(2);
eKm = &myMats->A[0];
eKg = &myMats->A[0];
myMats->rhsOnly = true;
break;
default:
this->Elasticity::setMode(mode);
return;
}
// We always need the force+displacement vectors in nonlinear simulations
iS = &myMats->b[0];
eS = &myMats->b[0];
eV = &myMats->b[1];
tracVal.clear();
}
bool NonlinearElasticityUL::evalInt (LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const
{
elmInt = myMats;
// Evaluate the deformation gradient, F, and the Green-Lagrange strains, E
SymmTensor E(nsd);
if (!this->kinematics(dNdX,E))
return false;
bool lHaveStrains = !E.isZero();
if (lHaveStrains)
{
// Invert the deformation gradient ==> Fi
Matrix Fi(F);
double J = Fi.inverse();
if (J == 0.0) return false;
// Scale with J=|F| since we are integrating over current configuration
detJW *= J;
if (eKm || iS)
{
// Push-forward the basis function gradients to current configuration
dNdx.multiply(dNdX,Fi); // dNdx = dNdX * F^-1
// Compute the small-deformation strain-displacement matrix B from dNdx
if (!this->formBmatrix(dNdx)) return false;
#if INT_DEBUG > 0
std::cout <<"NonlinearElasticityUL::dNdx ="<< dNdx;
std::cout <<"NonlinearElasticityUL::B ="<< Bmat;
#endif
}
}
else if (eKm || iS)
// Initial state, no deformation yet
if (!this->formBmatrix(dNdX)) return false;
// Evaluate the constitutive relation
SymmTensor sigma(nsd);
double U = 0.0;
if (eKm || eKg || iS)
if (!this->constitutive(Cmat,sigma,U,E,X, (eKg || iS) && lHaveStrains))
return false;
if (eKm)
{
// Integrate the material stiffness matrix
CB.multiply(Cmat,Bmat).multiply(detJW); // CB = C*B*|J|*w
eKm->multiply(Bmat,CB,true,false,true); // EK += B^T * CB
}
if (eKg && lHaveStrains)
// Integrate the geometric stiffness matrix
this->formKG(*eKg,dNdx,sigma,detJW);
if (eM)
// Integrate the mass matrix
this->formMassMatrix(*eM,N,X,detJW);
if (iS && lHaveStrains)
{
// Integrate the internal forces
sigma *= -detJW;
if (!Bmat.multiply(sigma,*iS,true,true)) // ES -= B^T*sigma
return false;
}
if (eS)
// Integrate the load vector due to gravitation and other body forces
this->formBodyForce(*eS,N,X,detJW);
return true;
}
bool NonlinearElasticityUL::evalBou (LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X, const Vec3& normal) const
{
elmInt = myMats;
if (!tracFld)
{
std::cerr <<" *** NonlinearElasticityUL::evalBou: No tractions."
<< std::endl;
return false;
}
else if (!eS)
{
std::cerr <<" *** NonlinearElasticityUL::evalBou: No load vector."
<< std::endl;
return false;
}
// Evaluate the surface traction
Vec3 T = (*tracFld)(X,normal);
// Store the traction value for vizualization
if (!T.isZero()) tracVal[X] = T;
unsigned short int i, j;
if (loadOp == 1)
{
// Compute the deformation gradient, F
if (!this->formDefGradient(dNdX)) return false;
// Check for with-rotated pressure load
if (tracFld->isNormalPressure())
{
// Compute its inverse and determinant, J
double J = F.inverse();
if (J == 0.0) return false;
// Pull-back the normal traction to the initial configuration.
// See equation (3.4.5) on page 102 in T. Belytschko et. al (2000):
// p*n*dS ==> J*p*n0*F^-1*dS0
Vec3 t = J*T; T = 0.0;
for (i = 1; i <= nsd; i++)
for (j = 1; j <= nsd; j++)
T[i-1] += t[j-1]*F(j,i);
}
else
// Scale with J=|F| since we are integrating over current configuration
detJW *= F.det();
}
// Integrate the force vector
Vector& ES = *eS;
for (size_t a = 1; a <= N.size(); a++)
for (i = 1; i <= nsd; i++)
ES(nsd*(a-1)+i) += T[i-1]*N(a)*detJW;
return true;
}
bool NonlinearElasticityUL::formDefGradient (const Matrix& dNdX) const
{
static SymmTensor dummy(0);
return this->kinematics(dNdX,dummy);
}
bool NonlinearElasticityUL::kinematics (const Matrix& dNdX, SymmTensor& E) const
{
if (!eV || eV->empty())
{
// Initial state, unit deformation gradient and zero strains
F.diag(1.0,nsd);
E.zero();
return true;
}
const size_t nenod = dNdX.rows();
const size_t nstrc = nsd*(nsd+1)/2;
if (eV->size() != nenod*nsd || dNdX.cols() < nsd)
{
std::cerr <<" *** NonlinearElasticityUL::kinematics: Invalid dimension,"
<<" dNdX("<< nenod <<","<< dNdX.cols() <<")"<< std::endl;
return false;
}
// Compute the deformation gradient, [F] = [I] + [dudX] = [I] + [dNdX]*[u],
// and the Green-Lagrange strains, E_ij = 0.5(F_ij+F_ji+F_ki*F_kj).
// Notice that the matrix multiplication method used here treats the
// element displacement vector, *eV, as a matrix whose number of columns
// equals the number of rows in the matrix dNdX.
if (!F.multiplyMat(*eV,dNdX)) // F = Grad{u} = eV*dNdX
return false;
unsigned short int i, j, k;
// Now form the Green-Lagrange strain tensor.
// Note that for the shear terms (i/=j) we actually compute 2*E_ij
// to be consistent with the engineering strain style constitutive matrix.
for (i = 1; i <= E.dim(); i++)
for (j = 1; j <= i; j++)
{
double Eij = F(i,j) + F(j,i);
for (k = 1; k <= nsd; k++)
Eij += F(k,i)*F(k,j);
E(i,j) = i == j ? 0.5*Eij : Eij;
}
// Add the unit tensor to F to form the deformation gradient
for (i = 1; i <= nsd; i++)
F(i,i) += 1.0;
#if INT_DEBUG > 0
std::cout <<"NonlinearElasticityUL::eV ="<< *eV;
std::cout <<"NonlinearElasticityUL::F ="<< F;
#endif
return true;
}
bool NonlinearElasticityUL::constitutive (Matrix& C, SymmTensor& sigma,
double& U, const SymmTensor& eps,
const Vec3& X,
char calcStress) const
{
double J = F.det();
if (J == 0.0)
{
std::cerr <<" *** NonlinearElasticityUL::constitutive: "
<<" Singular/zero deformation gradient"<< F;
return false;
}
size_t ndim = F.rows();
#ifdef USE_FTNMAT
if (mVER >= 0)
{
// Invoke the FORTRAN routine for constitutive relations.
// This also includes the Neo-Hookean hyperelastic material models.
size_t ncmp = ndim*(ndim+1)/2;
C.resize(ncmp,ncmp);
int ipsw = INT_DEBUG > 1 ? 9 : 0;
int ierr = 0;
if (ndim == 2)
cons2d_(mTYP,mVER,J,F.ptr(),pmat,U,sigma.ptr(),C.ptr(),ipsw,6,ierr);
else
cons3d_(mTYP,mVER,J,F.ptr(),pmat,U,sigma.ptr(),C.ptr(),ipsw,6,ierr);
if (calcStress == 2)
{
// Transform to 2nd Piola-Kirchhoff stresses,
// via pull-back to reference configuration
F.inverse();
Tensor Fi(ndim);
sigma.transform(Fi=F); // sigma = F^-1 * sigma * F^-t
sigma *= J;
}
#if INT_DEBUG > 0
if (calcStress)
std::cout <<"NonlinearElasticityUL::sigma =\n"<< sigma;
std::cout <<"NonlinearElasticityUL::C ="<< C;
#endif
return ierr == 0;
}
#endif
// We are only doing pure isotrophic linear elastic material in this method.
// Evaluate the constitutive matrix at this point, but first reset nsd to ndim
// since it is used to set the dimension on C.
unsigned short int tmp = nsd;
const_cast<NonlinearElasticityUL*>(this)->nsd = ndim;
if (!this->formCmatrix(C,X)) return false;
const_cast<NonlinearElasticityUL*>(this)->nsd = tmp;
if (calcStress)
// Evaluate the symmetric stress tensor
C.multiply(eps,sigma); // sigma = C*eps
if (calcStress == 1)
{
// Push-forward the stress tensor to current configuration
Tensor dUdX(ndim);
sigma.transform(dUdX=F); // sigma = F * sigma * F^t
sigma *= 1.0/J;
#if INT_DEBUG > 0
std::cout <<"NonlinearElasticityUL::sigma =\n"<< sigma;
#endif
}
else if (calcStress == 2)
return true;
// Push-forward the constitutive matrix to current configuration
size_t i, j, k, l, m = C.rows();
Matrix T(m,m), Ctmp;
for (i = 1; i <= ndim; i++)
{
for (j = 1; j <= ndim; j++)
T(i,j) = F(j,i)*F(j,i);
for (j = ndim+1, k = 1; j <= m; j++, k++)
T(i,j) = 0.5*(F(k,i)*F(k%3+1,i) + F(k%3+1,i)*F(k,i));
}
for (i = ndim+1, k = 1; i <= m; i++, k++)
{
for (j = 1; j <= ndim; j++)
T(i,j) = F(j,k)*F(j,k%3+1) + F(j,k%3+1)*F(j,k);
for (j = ndim+1, l = 1; j <= m; j++, l++)
T(i,j) = 0.5*(F(l,k)*F(l%3+1,k%3+1) + F(l%3+1,k)*F(l,k%3+1) +
F(l,k%3+1)*F(l%3+1,k) + F(l%3+1,k%3+1)*F(l,k));
}
// C = 1/J * T^t * C * T
Ctmp.multiply(C,T);
C.multiply(T,Ctmp,true);
C *= 1.0/J;
#if INT_DEBUG > 0
std::cout <<"NonlinearElasticityUL::C ="<< C;
#endif
return true;
}
NormBase* NonlinearElasticityUL::getNormIntegrand (TensorFunc*) const
{
return new ElasticityNormUL(*const_cast<NonlinearElasticityUL*>(this));
}
bool ElasticityNormUL::evalInt (LocalIntegral*& elmInt, double detJW,
const Vector&, const Matrix& dNdX,
const Vec3& X) const
{
ElmNorm& pnorm = ElasticityNorm::getElmNormBuffer(elmInt);
NonlinearElasticityUL* ulp = static_cast<NonlinearElasticityUL*>(&problem);
// Evaluate the deformation gradient, F, and the Green-Lagrange strains, E
SymmTensor E(dNdX.cols());
if (!ulp->kinematics(dNdX,E))
return false;
// Evaluate the 2nd Piola Kirchhoff stresses, S
Matrix C;
SymmTensor S(E.dim());
double U = 0.0;
if (!ulp->constitutive(C,S,U,E,X,2))
return false;
if (U == 0.0)
{
// Integrate the energy norm a(u^h,u^h) = Int_Omega0 (S:E) dV0
const RealArray& sigma = S;
const RealArray& epsil = E;
for (size_t i = 0; i < sigma.size(); i++)
U += sigma[i]*epsil[i];
}
pnorm[0] += U*detJW;
return true;
}

View File

@@ -0,0 +1,204 @@
// $Id: NonlinearElasticityUL.h,v 1.7 2011-02-08 09:06:02 kmo Exp $
//==============================================================================
//!
//! \file NonlinearElasticityUL.h
//!
//! \date Sep 21 2010
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Integrand implementations for nonlinear elasticity problems.
//!
//==============================================================================
#ifndef _NONLINEAR_ELASTICITY_UL_H
#define _NONLINEAR_ELASTICITY_UL_H
#include "Elasticity.h"
/*!
\brief Class representing the integrand of the nonlinear elasticity problem.
\details This class implements an Updated Lagrangian formulation. It inherits
most of the Elasticity methods, but reimplements the \a kinematics method
for calculating the deformation gradient and the associated Green-Lagrange
strain tensor, and the method \a constitutive for calculating the tangent
constitutive matrix and the associated stress tensor. The \a evalInt and
\a evalBou methods are also reimplemented to account for the updated geometry.
This class supports isotrophic linear-elastic and neo-Hookean hyperelastic
material laws. The material properties are assumed constant in space.
*/
class NonlinearElasticityUL : public Elasticity
{
public:
//! \brief The default constructor invokes the parent class constructor.
//! \param[in] n Number of spatial dimensions
//! \param[in] ps If \e true, assume plane stress in 2D
//! \param[in] mver Material version parameter
//! \param[in] lop Load option (0=on initial length, 1=on updated length)
NonlinearElasticityUL(unsigned short int n = 3,
bool ps = false, int mver = 0, char lop = 0);
//! \brief Empty destructor.
virtual ~NonlinearElasticityUL() {}
//! \brief Prints out problem definition to the given output stream.
virtual void print(std::ostream& os) const;
//! \brief Defines the material properties to use.
//! \param[in] Young Young's modulus
//! \param[in] Poiss Poisson's ratio
//! \param[in] Density Mass density
//!
//! \details If \a Poiss > 0.5, it is assumed that \a Young and \a Poiss
//! are the Lame's parameters instead (bulk and shear moduli).
virtual void setMaterial(double Young, double Poiss, double Density);
//! \brief Defines the solution mode before the element assembly is started.
//! \param[in] mode The solution mode to use
virtual void setMode(SIM::SolutionMode mode);
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] X Cartesian coordinates of current integration point
virtual bool evalInt(LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const;
//! \brief Evaluates the integrand at a boundary point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] X Cartesian coordinates of current integration point
//! \param[in] normal Boundary normal vector at current integration point
//!
//! \details This method is reimplemented in this class to account for
//! possibly with-rotated traction fields (non-conservative loads).
//! For uni-directional (conservative) loads, it is similar to the
//! \a LinearElasticity::evalBou method.
virtual bool evalBou(LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X, const Vec3& normal) const;
//! \brief Returns a pointer to an Integrand for solution norm evaluation.
//! \note The Integrand object is allocated dynamically and has to be deleted
//! manually when leaving the scope of the pointer variable receiving the
//! returned pointer value.
virtual NormBase* getNormIntegrand(TensorFunc* = 0) const;
//! \brief Calculates some kinematic quantities at current point.
//! \param[in] dNdX Basis function gradients at current point
//! \param[out] E Green-Lagrange strain tensor at current point
//!
//! \details The deformation gradient is established and stored in
//! the mutable class member \a F.
virtual bool kinematics(const Matrix& dNdX, SymmTensor& E) const;
//! \brief Calculates the deformation gradient at current point.
//! \param[in] dNdX Basis function gradients at current point
bool formDefGradient(const Matrix& dNdX) const;
//! \brief Evaluates the constitutive relation at current point.
//! \param[out] C Constitutive matrix at current point
//! \param[out] sigma Stress tensor at current point
//! \param[out] U Strain energy density
//! \param[in] eps Strain tensor at current point
//! \param[in] X Cartesian coordinates of current point
//! \param[in] calcStress Stress calculation option.
//! 0: No stress calculation, only output the C-matrix,
//! 1: Calculate Cauchy stresses,
//! 2: Calculate Second Piola Kirchhoff stresses.
virtual bool constitutive(Matrix& C, SymmTensor& sigma, double& U,
const SymmTensor& eps, const Vec3& X,
char calcStress = 1) const;
protected:
mutable Matrix F; //!< Deformation gradient
mutable Matrix dNdx; //!< Basis function gradients in current configuration
mutable Matrix CB; //!< Result of the matrix-matrix product C*B
private:
char loadOp; //!< Load option
int mTYP; //!< Material type
int mVER; //!< Material version
double pmat[2]; //!< Material properties
};
/*!
\brief Class representing the integrand of the elasticity energy norm.
\details This class reimplements the \a evalInt method to use the strain
energy density value returned by the nonlinear constitutive model.
It also forwards the mixed interface methods to the corresponding single-field
methods of the parent class.
*/
class ElasticityNormUL : public ElasticityNorm
{
public:
//! \brief The only constructor initializes its data members.
//! \param[in] p The linear elasticity problem to evaluate norms for
ElasticityNormUL(NonlinearElasticityUL& p) : ElasticityNorm(p,0) {}
//! \brief Empty destructor.
virtual ~ElasticityNormUL() {}
//! \brief Initializes current element for numerical integration (mixed).
//! \param[in] MNPC1 Nodal point correspondance for the displacement field
virtual bool initElement(const std::vector<int>& MNPC1,
const std::vector<int>&, size_t)
{
return this->ElasticityNorm::initElement(MNPC1);
}
//! \brief Initializes current element for boundary integration (mixed).
//! \param[in] MNPC1 Nodal point correspondance for the displacement field
virtual bool initElementBou(const std::vector<int>& MNPC1,
const std::vector<int>&, size_t)
{
return this->ElasticityNorm::initElement(MNPC1);
}
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] X Cartesian coordinates of current integration point
virtual bool evalInt(LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const;
//! \brief Evaluates the integrand at an interior point (mixed).
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N1 Basis function values for the displacement field
//! \param[in] dN1dX Basis function gradients for the displacement field
//! \param[in] X Cartesian coordinates of current integration point
virtual bool evalInt(LocalIntegral*& elmInt, double detJW,
const Vector& N1, const Vector&,
const Matrix& dN1dX, const Matrix&,
const Vec3& X) const
{
return this->evalInt(elmInt,detJW,N1,dN1dX,X);
}
//! \brief Evaluates the integrand at a boundary point (mixed).
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N1 Basis function values for the displacement field
//! \param[in] dN1dX Basis function gradients for the displacement field
//! \param[in] X Cartesian coordinates of current integration point
//! \param[in] normal Boundary normal vector at current integration point
virtual bool evalBou(LocalIntegral*& elmInt, double detJW,
const Vector& N1, const Vector&,
const Matrix& dN1dX, const Matrix&,
const Vec3& X, const Vec3& normal) const
{
return this->ElasticityNorm::evalBou(elmInt,detJW,N1,dN1dX,X,normal);
}
};
#endif

View File

@@ -0,0 +1,595 @@
// $Id: NonlinearElasticityULMX.C,v 1.7 2011-02-08 09:06:02 kmo Exp $
//==============================================================================
//!
//! \file NonlinearElasticityULMX.C
//!
//! \date Dec 14 2010
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Integrand implementations for nonlinear elasticity mixed problems.
//!
//==============================================================================
#include "NonlinearElasticityULMX.h"
#include "ElmMats.h"
#include "ElmNorm.h"
#include "Utilities.h"
#include "Vec3Oper.h"
#ifdef USE_FTNMAT
extern "C" {
//! \brief Project the constitutive matrix for mixed models.
void pcst3d_(const double* C, double* D,
const int& ipsw, const int& iwr);
//! \brief Compute finite deformation mixed tangent material matrix.
void mdma3d_(const double& P_bar, const double& P_mix,
const double* Sig, double* D,
const int& ipsw, const int& iwr);
//! \brief Accumulates material stiffness contributions for mixed 2D problems.
void acckmx2d_(const int& nEN, const double* dNdx, const double* dNdxBar,
const double* D, double* eKt);
//! \brief Accumulates material stiffness contributions for mixed 3D problems.
void acckmx3d_(const int& nEN, const double* dNdx, const double* dNdxBar,
const double* D, double* eKt);
}
#ifndef INT_DEBUG
#define INT_DEBUG 0
#endif
#endif
NonlinearElasticityULMX::NonlinearElasticityULMX (unsigned short int n,
int mver, int pp)
: NonlinearElasticityUL(n, false, mver < 0 ? 0 : mver)
{
p = pp;
iP = 0;
// Both the current and previous solutions are needed
primsol.resize(2);
}
void NonlinearElasticityULMX::print (std::ostream& os) const
{
this->NonlinearElasticityUL::print(os);
std::cout <<"NonlinearElasticityULMX: Mixed formulation,"
<<" discontinuous pressure, p="<< p << std::endl;
}
void NonlinearElasticityULMX::setMode (SIM::SolutionMode mode)
{
if (!myMats) return;
myMats->rhsOnly = false;
Hh = eM = eKm = eKg = 0;
iS = eS = eV = 0;
switch (mode)
{
case SIM::STATIC:
myMats->resize(2,3);
eKm = &myMats->A[0];
eKg = &myMats->A[0];
Hh = &myMats->A[1];
break;
case SIM::DYNAMIC:
myMats->resize(3,3);
eKm = &myMats->A[0];
eKg = &myMats->A[0];
eM = &myMats->A[1];
Hh = &myMats->A[2];
break;
case SIM::RHS_ONLY:
if (myMats->A.size() < 2)
myMats->resize(2,3);
else
myMats->b.resize(3);
eKm = &myMats->A[0];
eKg = &myMats->A[0];
Hh = &myMats->A.back();
myMats->rhsOnly = true;
break;
case SIM::RECOVERY:
if (myMats->A.size() < 1)
myMats->resize(1,3);
else
myMats->b.resize(3);
Hh = &myMats->A.back();
eV = &myMats->b[1];
myMats->rhsOnly = true;
return;
default:
this->Elasticity::setMode(mode);
return;
}
// We always need the force+displacement vectors in nonlinear simulations
iS = &myMats->b[0];
eS = &myMats->b[0];
eV = &myMats->b[1];
tracVal.clear();
}
bool NonlinearElasticityULMX::initElement (const std::vector<int>& MNPC,
const Vec3& Xc, size_t nPt)
{
#if INT_DEBUG > 0
std::cout <<"\n\n *** Entering NonlinearElasticityULMX::initElement";
std::cout <<", Xc = "<< Xc << std::endl;
#endif
iP = 0;
X0 = Xc;
myData.resize(nPt);
size_t nPm = utl::Pascal(p,nsd);
if (Hh) Hh->resize(nPm,nPm,true);
// Extract the element level displacements of previous increment
int ierr = 0;
if (!primsol[1].empty())
if (ierr = utl::gather(MNPC,nsd,primsol[1],myMats->b[2]))
std::cerr <<" *** NonlinearElasticityULMX::initElement: Detected "
<< ierr <<" node numbers out of range."<< std::endl;
// The other element matrices are initialized by the parent class method
return this->NonlinearElasticityUL::initElement(MNPC) && ierr == 0;
}
bool NonlinearElasticityULMX::evalInt (LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const
{
if (myMats->b.size() < 3) return false;
ItgPtData& ptData = myData[iP++];
unsigned short int i, j;
#if INT_DEBUG > 0
std::cout <<"NonlinearElasticityUL::dNdX ="<< dNdX;
#endif
// Evaluate the deformation gradient, Fp, at previous configuration
const_cast<NonlinearElasticityULMX*>(this)->eV = &myMats->b[2];
if (this->formDefGradient(dNdX))
for (i = 1; i <= nsd; i++)
for (j = 1; j <= nsd; j++)
ptData.Fp(i,j) = F(i,j);
else
return false;
// Evaluate the deformation gradient, F, at current configuration
const_cast<NonlinearElasticityULMX*>(this)->eV = &myMats->b[1];
if (this->formDefGradient(dNdX))
for (i = 1; i <= nsd; i++)
for (j = 1; j <= nsd; j++)
ptData.F(i,j) = F(i,j);
else
return false;
if (nsd == 2) // In 2D we always assume plane strain so set F(3,3)=1
ptData.F(3,3) = ptData.Fp(3,3) = 1.0;
// Invert the deformation gradient ==> Fi
double J = F.inverse();
if (J == 0.0) return false;
// Push-forward the basis function gradients to current configuration
ptData.dNdx.multiply(dNdX,F); // dNdx = dNdX * F^-1
ptData.X.assign(X);
ptData.detJW = detJW;
// Evaluate the pressure modes (generalized coordinates)
Vec3 Xg = X-X0;
if (nsd == 3)
utl::Pascal(p,Xg.x,Xg.y,Xg.z,ptData.Phi);
else if (nsd == 2)
utl::Pascal(p,Xg.x,Xg.y,ptData.Phi);
else
return false;
if (Hh)
// Integrate the Hh-matrix
Hh->outer_product(ptData.Phi,ptData.Phi*detJW,true);
if (eM)
// Integrate the mass matrix
this->formMassMatrix(*eM,N,X,J*detJW);
if (eS)
// Integrate the load vector due to gravitation and other body forces
this->formBodyForce(*eS,N,X,J*detJW);
return this->getIntegralResult(elmInt);
}
bool NonlinearElasticityULMX::finalizeElement (LocalIntegral*& elmInt)
{
if (!iS && !eKm && !eKg) return true;
if (!Hh) return false;
#if INT_DEBUG > 0
std::cout <<"\n\n *** Entering NonlinearElasticityULMX::finalizeElement\n";
for (iP = 1; iP <= myData.size(); iP++)
{
const ItgPtData& pt = myData[iP-1];
std::cout <<"\n X" << iP <<" = "<< pt.X;
std::cout <<"\n detJ"<< iP <<"W = "<< pt.detJW;
std::cout <<"\n F" << iP <<" =\n"<< pt.F;
std::cout <<" dNdx"<< iP <<" ="<< pt.dNdx;
std::cout <<" Phi" << iP <<" ="<< pt.Phi;
}
std::cout <<"\n H ="<< *Hh << std::endl;
#endif
// 1. Eliminate the internal pressure DOFs by static condensation.
size_t i, j;
size_t nPM = Hh->rows();
size_t nEN = myData.front().dNdx.rows();
size_t nGP = myData.size();
// Invert the H-matrix
if (!utl::invert(*Hh)) return false;
const Matrix& Hi = *Hh;
Matrix Theta(2,nGP), dNdxBar[nGP];
if (nPM == 1)
{
// Simplified calculation when the only pressure mode is Phi=1.0
double h1 = 0.0;
double h2 = 0.0;
dNdxBar->resize(nEN,nsd,true);
for (iP = 0; iP < nGP; iP++)
{
const ItgPtData& pt = myData[iP];
h1 += pt.F.det() * pt.detJW;
h2 += pt.Fp.det()* pt.detJW;
dNdxBar->add(pt.dNdx,pt.F.det()*pt.detJW);
}
dNdxBar->multiply(1.0/h1);
// All gauss point values are identical
Theta(1,1) = h1 * Hi(1,1);
Theta(2,1) = h2 * Hi(1,1);
for (iP = 2; iP <= nGP; iP++)
{
Theta(1,iP) = Theta(1,1);
Theta(2,iP) = Theta(2,1);
dNdxBar[iP-1] = dNdxBar[0];
}
}
else // higher order elements with nPM > 1
{
Matrix Gg[nPM], Hg[nPM];
Vector J1(nPM), J2(nPM), Hj1, Hj2;
for (i = 0; i < nPM; i++)
{
Gg[i].resize(nEN,nsd,true);
Hg[i].resize(nEN,nsd,true);
}
// Integrate the Ji- and G-matrices
for (iP = 0; iP < nGP; iP++)
{
const ItgPtData& pt = myData[iP];
for (i = 0; i < nPM; i++)
{
double h1 = pt.Phi[i] * pt.F.det() * pt.detJW;
double h2 = pt.Phi[i] * pt.Fp.det()* pt.detJW;
J1[i] += h1;
J2[i] += h2;
Gg[i].add(pt.dNdx,h1);
}
}
#if INT_DEBUG > 0
std::cout <<"NonlinearElasticityULMX::J1 ="<< J1;
std::cout <<"NonlinearElasticityULMX::J2 ="<< J2;
for (i = 0; i < nPM; i++)
std::cout <<"NonlinearElasticityULMX::G"<< i+1 <<" ="<< Gg[i];
#endif
// Calculate Hji = Hi*Ji
if (!Hi.multiply(J1,Hj1)) return false;
if (!Hi.multiply(J2,Hj2)) return false;
// Calculate Hg = Hi*Gg
for (i = 1; i <= nPM; i++)
for (j = 1; j <= nPM; j++)
Hg[i-1].add(Gg[j-1],Hi(i,j));
#if INT_DEBUG > 0
std::cout <<"NonlinearElasticityULMX::Hj1 ="<< Hj1;
std::cout <<"NonlinearElasticityULMX::Hj2 ="<< Hj2;
for (i = 0; i < nPM; i++)
std::cout <<"NonlinearElasticityULMX::Hg"<< i+1 <<" ="<< Hg[i];
#endif
for (iP = 0; iP < nGP; iP++)
{
i = iP+1;
ItgPtData& pt = myData[iP];
// Calculate Theta = Hj*Phi
Theta(1,i) = Hj1.dot(pt.Phi);
Theta(2,i) = Hj2.dot(pt.Phi);
// Calculate dNdxBar = Hg*Phi * 1/Theta
dNdxBar[iP].resize(nEN,nsd,true);
for (j = 0; j < nPM; j++)
dNdxBar[iP].add(Hg[j],pt.Phi[j]/Theta(1,i));
}
}
// Modify the deformation gradients
Vector detF(nGP);
for (iP = 1; iP <= nGP; iP++)
{
ItgPtData& pt = myData[iP-1];
detF(iP) = pt.F.det();
pt.F *= pow(fabs(Theta(1,iP)/detF(iP)),1.0/3.0);
pt.Fp *= pow(fabs(Theta(2,iP)/pt.Fp.det()),1.0/3.0);
}
#if INT_DEBUG > 0
std::cout <<"NonlinearElasticityULMX::detF ="<< detF;
std::cout <<"NonlinearElasticityULMX::Theta ="<< Theta;
for (iP = 0; iP < nGP; iP++)
std::cout <<"NonlinearElasticityULMX::dNdxBar"<< iP+1 <<" ="<< dNdxBar[iP];
#endif
// 2. Evaluate the constitutive relation and calculate mixed pressure state.
Matrix D[nGP];
Vector Sigm(nPM), Hsig;
std::vector<SymmTensor> Sig;
Sig.reserve(nGP);
F.resize(3,3);
bool lHaveStress = false;
for (iP = 0; iP < nGP; iP++)
{
#if INT_DEBUG > 0
std::cout <<"\n iGP = "<< iP+1 << std::endl;
#endif
const ItgPtData& pt = myData[iP];
// Evaluate the constitutive relation
double U = 0.0;
F.fill(pt.F.ptr());
Sig.push_back(SymmTensor(3));
if (!this->constitutive(Cmat,Sig.back(),U,Sig.back(),pt.X))
return false;
#ifdef USE_FTNMAT
// Project the constitutive matrix for the mixed model
D[iP].resize(7,7);
int ipsw = INT_DEBUG > 1 ? 9 : 0;
pcst3d_(Cmat.ptr(),D[iP].ptr(),ipsw,6);
#else
std::cerr <<" *** NonlinearElasticityULMX::finalizeElement: "
<<" Does not work when compiled without USE_FTNMAT"<< std::endl;
return false;
#endif
if (U > 0.0)
{
// Integrate mean pressure
lHaveStress = true;
double Mpress = Sig.back().trace()*pt.detJW/3.0;
for (i = 1; i <= nPM; i++)
Sigm(i) += Mpress * pt.Phi(i);
}
}
// Divide pressure by reference volume; Hsig = Hi*Sigm
if (lHaveStress)
if (!Hh->multiply(Sigm,Hsig)) return false;
// 3. Evaluate tangent matrices.
double Press = 0.0, Bpres = 0.0;
SymmTensor Sigma(nsd);
for (iP = 0; iP < nGP; iP++)
{
const ItgPtData& pt = myData[iP];
double dFmix = Theta(1,iP+1);
if (lHaveStress)
{
Press = Hsig.dot(pt.Phi) * detF[iP]/dFmix;
Bpres = Sig[iP].trace()/3.0;
Sigma = Sig[iP] + (Press-Bpres);
#if INT_DEBUG > 0
std::cout <<"NonlinearElasticityULMX::Sigma"<< iP+1 <<"\n"<< Sigma;
#endif
if (eKg)
// Integrate the geometric stiffness matrix
this->formKG(*eKg,pt.dNdx,Sigma,dFmix*pt.detJW);
}
if (eKm)
{
#ifdef USE_FTNMAT
int ipsw = INT_DEBUG > 1 ? 9 : 0;
mdma3d_(Bpres,Press,Sig[iP].ptr(),D[iP].ptr(),ipsw,6);
// Integrate the material stiffness matrix
D[iP] *= dFmix*pt.detJW;
if (nsd == 2)
acckmx2d_(nEN,pt.dNdx.ptr(),dNdxBar[iP].ptr(),D[iP].ptr(),eKm->ptr());
else
acckmx3d_(nEN,pt.dNdx.ptr(),dNdxBar[iP].ptr(),D[iP].ptr(),eKm->ptr());
#endif
}
if (iS && lHaveStress)
{
// Compute the small-deformation strain-displacement matrix B from dNdx
if (!this->formBmatrix(pt.dNdx)) return false;
// Integrate the internal forces
Sigma *= -dFmix*pt.detJW;
if (!Bmat.multiply(Sigma,*iS,true,true)) // ES -= B^T*sigma
return false;
}
}
#if INT_DEBUG > 0
std::cout <<"\n *** Leaving NonlinearElasticityULMX::finalizeElement";
if (eKm) std::cout <<", eKm ="<< *eKm;
if (iS) std::cout <<" iS ="<< *iS;
std::cout << std::endl;
#endif
return true;
}
NormBase* NonlinearElasticityULMX::getNormIntegrand (TensorFunc*) const
{
return new ElasticityNormULMX(*const_cast<NonlinearElasticityULMX*>(this));
}
bool ElasticityNormULMX::initElement (const std::vector<int>& MNPC,
const Vec3& Xc, size_t nPt)
{
NonlinearElasticityULMX& elp = static_cast<NonlinearElasticityULMX&>(problem);
elp.iP = 0;
elp.X0 = Xc;
elp.myData.resize(nPt);
size_t nPm = utl::Pascal(elp.p,elp.nsd);
if (elp.Hh) elp.Hh->resize(nPm,nPm,true);
// The other element matrices are initialized by the parent class method
return elp.NonlinearElasticityUL::initElement(MNPC);
}
bool ElasticityNormULMX::evalInt (LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const
{
NonlinearElasticityULMX& elp = static_cast<NonlinearElasticityULMX&>(problem);
return elp.evalInt(elmInt,detJW,N,dNdX,X);
}
bool ElasticityNormULMX::finalizeElement (LocalIntegral*& elmInt)
{
NonlinearElasticityULMX& elp = static_cast<NonlinearElasticityULMX&>(problem);
if (!elp.Hh) return false;
size_t iP;
#if INT_DEBUG > 0
std::cout <<"\n\n *** Entering ElasticityNormULMX::finalizeElement\n";
for (iP = 1; iP <= elp.myData.size(); iP++)
{
const NonlinearElasticityULMX::ItgPtData& pt = elp.myData[iP-1];
std::cout <<"\n X" << iP <<" = "<< pt.X;
std::cout <<"\n detJ"<< iP <<"W = "<< pt.detJW;
std::cout <<"\n F" << iP <<" =\n"<< pt.F;
std::cout <<" dNdx"<< iP <<" ="<< pt.dNdx;
std::cout <<" Phi" << iP <<" ="<< pt.Phi;
}
std::cout <<"\n H ="<< *(elp.Hh) << std::endl;
#endif
// 1. Eliminate the internal pressure DOFs by static condensation.
size_t nPM = elp.Hh->rows();
size_t nGP = elp.myData.size();
// Invert the H-matrix
Matrix& Hi = *elp.Hh;
if (!utl::invert(Hi)) return false;
Vector Theta(nGP);
if (nPM == 1)
{
// Simplified calculation when the only pressure mode is Phi=1.0
double h1 = 0.0;
for (iP = 0; iP < nGP; iP++)
h1 += elp.myData[iP].F.det() * elp.myData[iP].detJW;
// All gauss point values are identical
Theta.front() = h1 * Hi(1,1);
for (iP = 1; iP < nGP; iP++)
Theta[iP] = Theta.front();
}
else // higher order elements with nPM > 1
{
Vector Ji(nPM), Hj;
// Integrate the Ji-matrix
for (iP = 0; iP < nGP; iP++)
{
const NonlinearElasticityULMX::ItgPtData& pt = elp.myData[iP];
for (size_t i = 0; i < nPM; i++)
Ji[i] += pt.Phi[i] * pt.F.det() * pt.detJW;
}
// Calculate Hj = Hi*Ji
if (!Hi.multiply(Ji,Hj)) return false;
#if INT_DEBUG > 0
std::cout <<"ElasticityNormULMX::Ji ="<< Ji;
std::cout <<"ElasticityNormULMX::Hj ="<< Hj;
#endif
// Calculate Theta = Hj*Phi
for (iP = 0; iP < nGP; iP++)
Theta[iP] = Hj.dot(elp.myData[iP].Phi);
}
#if INT_DEBUG > 0
std::cout <<"ElasticityNormULMX::Theta ="<< Theta;
#endif
// 2. Evaluate the constitutive relation and integrate energy.
ElmNorm& pnorm = ElasticityNorm::getElmNormBuffer(elmInt);
Matrix C;
SymmTensor Sig(3);
elp.F.resize(3,3);
for (iP = 0; iP < nGP; iP++)
{
#if INT_DEBUG > 0
std::cout <<"\n iGP = "<< iP+1 << std::endl;
#endif
NonlinearElasticityULMX::ItgPtData& pt = elp.myData[iP];
// Modify the deformation gradient
pt.F *= pow(fabs(Theta[iP]/pt.F.det()),1.0/3.0);
// Compute the strain energy density
double U = 0.0;
elp.F.fill(pt.F.ptr());
if (!elp.constitutive(C,Sig,U,Sig,pt.X,false))
return false;
// Integrate strain energy
pnorm[0] += U*pt.detJW;
}
return true;
}

View File

@@ -0,0 +1,156 @@
// $Id: NonlinearElasticityULMX.h,v 1.3 2011-01-08 16:29:10 kmo Exp $
//==============================================================================
//!
//! \file NonlinearElasticityULMX.h
//!
//! \date Dec 14 2010
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Integrand implementations for nonlinear elasticity mixed problems.
//!
//==============================================================================
#ifndef _NONLINEAR_ELASTICITY_UL_MX_H
#define _NONLINEAR_ELASTICITY_UL_MX_H
#include "NonlinearElasticityUL.h"
#include "Tensor.h"
/*!
\brief Class representing the integrand of the nonlinear elasticity problem.
\details This class implements a mixed Updated Lagrangian formulation with
internal pressure and volumetric change modes.
*/
class NonlinearElasticityULMX : public NonlinearElasticityUL
{
//! \brief A struct with integration point data needed by \a finalizeElement.
struct ItgPtData
{
Tensor F; //!< Deformation gradient at current step/iteration
Tensor Fp; //!< Deformation gradient at previous (converged) step
Matrix dNdx; //!< Basis function gradients at current configuration
Vector Phi; //!< Internal modes for the pressure/volumetric-change fields
Vec4 X; //!< Cartesian coordinates of current integration point
double detJW; //!< Jacobian determinant times integration point weight
//! \brief Default constructor.
ItgPtData() : F(3), Fp(3) { detJW = 0.0; }
};
public:
//! \brief The default constructor invokes the parent class constructor.
//! \param[in] n Number of spatial dimensions
//! \param[in] mver Material version parameter
//! \param[in] pp Polynomial order of the pressure/volumetric-change field
NonlinearElasticityULMX(unsigned short int n = 3, int mver = 0, int pp = 1);
//! \brief Empty destructor.
virtual ~NonlinearElasticityULMX() {}
//! \brief Prints out problem definition to the given output stream.
virtual void print(std::ostream& os) const;
//! \brief Defines the solution mode before the element assembly is started.
//! \param[in] mode The solution mode to use
virtual void setMode(SIM::SolutionMode mode);
//! \brief Initializes current element for numerical integration.
//! \param[in] MNPC Matrix of nodal point correspondance for current element
//! \param[in] Xc Cartesian coordinates of the element center
//! \param[in] nPt Number of integration points in this element
virtual bool initElement(const std::vector<int>& MNPC,
const Vec3& Xc, size_t nPt);
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] X Cartesian coordinates of current integration point
//!
//! \details This method mainly stores the integration point quantities
//! depending on the basis function values in the internal member \a myData.
//! The actual numerical integration of the tangent stiffness and internal
//! forces is then performed by the \a finalizeElement method.
virtual bool evalInt(LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const;
//! \brief Finalizes the element matrices after the numerical integration.
//! \details This method is used to implement the multiple integration point
//! loops within an element, required by the mixed formulation with internal
//! modes. All data needed during the integration has been stored internally
//! in the \a myData member by the \a evalInt method.
//! \param elmInt The local integral object to receive the contributions
virtual bool finalizeElement(LocalIntegral*&);
//! \brief Returns a pointer to an Integrand for solution norm evaluation.
//! \note The Integrand object is allocated dynamically and has to be deleted
//! manually when leaving the scope of the pointer variable receiving the
//! returned pointer value.
virtual NormBase* getNormIntegrand(TensorFunc* = 0) const;
//! \brief Returns which integrand to be used.
virtual int getIntegrandType() const { return 4; }
private:
int p; //!< Polynomial order of the internal pressure field
Vec3 X0; //!< Cartesian coordinates of the element center
Matrix* Hh; //!< Pointer to element Hh-matrix associated with pressure modes
mutable size_t iP; //!< Local integration point counter
mutable std::vector<ItgPtData> myData; //!< Local integration point data
friend class ElasticityNormULMX;
};
/*!
\brief Class representing the integrand of the mixed elasticity energy norm.
*/
class ElasticityNormULMX : public ElasticityNormUL
{
public:
//! \brief The only constructor initializes its data members.
//! \param[in] p The linear elasticity problem to evaluate norms for
ElasticityNormULMX(NonlinearElasticityULMX& p) : ElasticityNormUL(p) {}
//! \brief Empty destructor.
virtual ~ElasticityNormULMX() {}
//! \brief Initializes current element for numerical integration.
//! \param[in] MNPC Matrix of nodal point correspondance for current element
//! \param[in] Xc Cartesian coordinates of the element center
//! \param[in] nPt Number of integration points in this element
virtual bool initElement(const std::vector<int>& MNPC,
const Vec3& Xc, size_t nPt);
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] X Cartesian coordinates of current integration point
//!
//! \details This method just forwards to the \a evalInt method of the
//! corresponding NonlinearElasticityULMX object, for which this class is
//! declared as friend such that it can access the data members. The actual
//! norm integration us then performed by the \a finalizeElement method.
virtual bool evalInt(LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const;
//! \brief Finalizes the element norms after the numerical integration.
//! \details This method is used to implement the multiple integration point
//! loops within an element required by the mixed formulation. It performs
//! a subset of the tasks done by NonlinearElasticityULMX::finalizeElement.
//! \param elmInt The local integral object to receive the contributions
virtual bool finalizeElement(LocalIntegral*& elmInt);
//! \brief Returns which integrand to be used.
virtual int getIntegrandType() const { return 4; }
};
#endif

View File

@@ -0,0 +1,435 @@
// $Id: NonlinearElasticityULMixed.C,v 1.4 2011-02-08 09:06:02 kmo Exp $
//==============================================================================
//!
//! \file NonlinearElasticityULMixed.C
//!
//! \date Dec 22 2010
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Integrand implementations for nonlinear elasticity mixed problems.
//!
//==============================================================================
#include "NonlinearElasticityULMixed.h"
#include "ElmMats.h"
#include "Utilities.h"
#include "Tensor.h"
#ifdef USE_FTNMAT
extern "C" {
//! \brief Project the constitutive matrix for mixed models.
void pcst3d_(const double* C, double* D,
const int& ipsw, const int& iwr);
//! \brief Compute finite deformation mixed tangent material matrix.
void mdma3d_(const double& P_bar, const double& P_mix,
const double* Sig, double* D,
const int& ipsw, const int& iwr);
//! \brief Accumulates material stiffness contributions for 2D problems.
void acckm2d_(const int& nEN, const double* dNdx,
const double* D, double* eKt);
//! \brief Accumulates material stiffness contributions for 3D problems.
void acckm3d_(const int& nEN, const double* dNdx,
const double* D, double* eKt);
}
#ifndef INT_DEBUG
#define INT_DEBUG 0
#endif
#endif
//! \brief Enum for element-level solution vectors
enum DisplacementVectors {
U = 0, // displacement
T = 1, // volumetric change (theta)
P = 2 // pressure
};
//! \brief Enum for element-level right-hand-side vectors
enum ResidualVectors {
Ru = 3,
Rt = 4,
Rp = 5,
Rres = 6,
NVEC = 7
};
//! \brief Enum for element-level left-hand-side matrices
enum TangentMatrices {
Kuu = 0,
Kut = 1,
Kup = 2,
Ktt = 3,
Ktp = 4,
Ktan = 5,
NMAT = 6
};
NonlinearElasticityULMixed::MixedElmMats::MixedElmMats ()
{
this->resize(NMAT,NVEC);
}
/*!
Tangent matrix for the mixed formulation:
\f[{\bf K}_{\rm T} = \left[\begin{array}{ccc}
{\bf K}_{uu} & {\bf K}_{u\theta} & {\bf K}_{up} \\
{\bf K}_{u\theta}^T & {\bf K}_{\theta\theta} & {\bf K}_{\theta p} \\
{\bf K}_{up}^T & {\bf K}_{\theta p}^T & {\bf 0}
\end{array}\right]\f]
where
-# \f${\bf K}_{uu}\f$ :
tangent stiffness matrix of the displacement formulation
-# \f${\bf K}_{\theta\theta}\f$ :
tangent stiffness matrix related to volumetric change
-# \f${\bf K}_{u\theta}\f$ :
coupling stiffness between displacements and volumetric change
-# \f${\bf K}_{up}\f$ :
coupling stiffness between displacements and pressure
-# \f${\bf K}_{\theta p}\f$ :
coupling stiffness between volumetric change and pressure
*/
const Matrix& NonlinearElasticityULMixed::MixedElmMats::getNewtonMatrix () const
{
#if INT_DEBUG > 0
std::cout <<"\nNonlinearElasticityULMixed::MixedElmMats::getNewtonMatrix:"
<<"\nKuu ="<< A[Kuu]
<<"\nKut ="<< A[Kut]
<<"\nKup ="<< A[Kup]
<<"\nKtt ="<< A[Ktt]
<<"\nKtp ="<< A[Ktp];
#endif
Matrix& N = const_cast<Matrix&>(A.back());
size_t i, j;
size_t n = A[Kuu].rows();
size_t m = A[Ktt].rows();
for (i = 1; i <= n; i++)
{
for (j = 1; j <= n; j++)
N(i,j) = A[Kuu](i,j);
for (j = 1; j <= m; j++)
{
size_t k = n+2*j-1;
size_t l = n+2*j;
N(i,k) = A[Kut](i,j);
N(k,i) = A[Kut](i,j);
N(i,l) = A[Kup](i,j);
N(l,i) = A[Kup](i,j);
}
}
for (i = 1; i <= m; i++)
for (j = 1; j <= m; j++)
{
size_t ki = n+2*i-1;
size_t kj = n+2*j-1;
size_t lj = n+2*j;
N(ki,kj) = A[Ktt](i,j);
N(ki,lj) = A[Ktp](i,j);
N(lj,ki) = A[Ktp](i,j);
}
#if INT_DEBUG > 0
std::cout <<"\nNewton matrix ="<< A.back();
#endif
return A.back();
}
/*!
Right-hand-side vector for the mixed formulation:
\f[{\bf R} = \left\{\begin{array}{c}
{\bf R}_u \\ {\bf R}_\theta \\ {\bf R}_p
\end{array}\right\}\f]
where
-# \f${\bf R}_{u}\f$ :
residual forces of the displacement DOFs
-# \f${\bf R}_\theta\f$ :
residual forces of the volumetric change DOFs
-# \f${\bf R}_p\f$ :
residual forces of the pressure DOFs
*/
const Vector& NonlinearElasticityULMixed::MixedElmMats::getRHSVector () const
{
#if INT_DEBUG > 0
std::cout <<"\nNonlinearElasticityULMixed::MixedElmMats::getRHSVector:"
<<"\nRu ="<< b[Ru];
if (withLHS)
std::cout <<"\nRt ="<< b[Rt] <<"\nRp ="<< b[Rp];
#endif
Vector& R = const_cast<Vector&>(b.back());
size_t i;
size_t n = b[Ru].size();
size_t m = b[Rt].size();
for (i = 1; i <= n; i++)
R(i) = b[Ru](i);
if (withLHS)
for (i = 1; i <= m; i++)
{
R(n+2*i-1) = b[Rt](i);
R(n+2*i ) = b[Rp](i);
}
else
std::fill(R.begin()+n,R.end(),0.0);
#if INT_DEBUG > 0
std::cout <<"\nRHS ="<< b.back();
#endif
return b.back();
}
NonlinearElasticityULMixed::NonlinearElasticityULMixed (unsigned short int n,
int mver)
: NonlinearElasticityUL(n,false,mver)
{
Fbar.resize(3,3);
Dmat.resize(7,7);
if (myMats) delete myMats;
myMats = new MixedElmMats();
eKm = &myMats->A[Kuu];
eKg = &myMats->A[Kuu];
iS = &myMats->b[Ru];
eS = &myMats->b[Ru];
eV = &myMats->b[U];
}
void NonlinearElasticityULMixed::print (std::ostream& os) const
{
this->NonlinearElasticityUL::print(os);
std::cout <<"NonlinearElasticityULMixed: "
<<"Continuous volumetric change and pressure fields"<< std::endl;
}
void NonlinearElasticityULMixed::setMode (SIM::SolutionMode mode)
{
switch (mode)
{
case SIM::INIT:
case SIM::STATIC:
tracVal.clear();
myMats->rhsOnly = false;
break;
case SIM::RHS_ONLY:
tracVal.clear();
case SIM::RECOVERY:
myMats->rhsOnly = true;
break;
default:
std::cerr <<"\n *** NonlinearElasticityULMixed::setMode: Invalid mode "
<< mode << std::endl;
}
}
bool NonlinearElasticityULMixed::initElement (const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2,
size_t n1)
{
const size_t nen1 = MNPC1.size();
const size_t nen2 = MNPC2.size();
const size_t nedof = nsd*nen1 + 2*nen2;
myMats->A[Kut].resize(nsd*nen1,nen2,true);
myMats->A[Kup].resize(nsd*nen1,nen2,true);
myMats->A[Ktt].resize(nen2,nen2,true);
myMats->A[Ktp].resize(nen2,nen2,true);
myMats->A.back().resize(nedof,nedof,true);
myMats->b[Rt].resize(nen2,true);
myMats->b[Rp].resize(nen2,true);
myMats->b.back().resize(nedof,true);
// Extract the element level volumetric change and pressure vectors
int ierr = 0;
if (!primsol.front().empty())
if (ierr = (utl::gather(MNPC2,0,2,primsol.front(),myMats->b[T],nsd*n1,n1) +
utl::gather(MNPC2,1,2,primsol.front(),myMats->b[P],nsd*n1,n1)))
std::cerr <<" *** NonlinearElasticityULMixed::initElement: Detected "
<< ierr/2 <<" node numbers out of range."<< std::endl;
// The other element matrices are initialized by the parent class method
return this->NonlinearElasticityUL::initElement(MNPC1) && ierr == 0;
}
bool NonlinearElasticityULMixed::initElementBou (const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2,
size_t)
{
const size_t nen1 = MNPC1.size();
const size_t nen2 = MNPC2.size();
const size_t nedof = nsd*nen1 + 2*nen2;
myMats->b[Ru].resize(nsd*nen1,true);
myMats->b.back().resize(nedof,true);
// Extract the element level displacement vector
int ierr = 0;
if (!primsol.front().empty())
if (ierr = utl::gather(MNPC1,nsd,primsol.front(),myMats->b[U]))
std::cerr <<" *** NonlinearElasticityULMixed::initElementBou: Detected "
<< ierr <<" node numbers out of range."<< std::endl;
myMats->withLHS = false;
return ierr == 0;
}
bool NonlinearElasticityULMixed::evalInt (LocalIntegral*& elmInt, double detJW,
const Vector& N1,
const Vector& N2,
const Matrix& dNdX1,
const Matrix& dNdX2,
const Vec3& X) const
{
// Evaluate the deformation gradient, F, and the Green-Lagrange strains, E
SymmTensor E(nsd);
if (!this->kinematics(dNdX1,E))
return false;
bool lHaveStrains = !E.isZero();
// Evaluate the volumetric change and pressure fields
double Theta = N2.dot(myMats->b[T]) + 1.0;
double Press = N2.dot(myMats->b[P]);
#if INT_DEBUG > 0
std::cout <<"NonlinearElasticityULMixed::b_theta ="<< myMats->b[T];
std::cout <<"NonlinearElasticityULMixed::b_p ="<< myMats->b[P];
std::cout <<"Theta = "<< Theta <<" Press = "<< Press << std::endl;
#endif
// Compute the mixed integration point volume
double dVol = Theta*detJW;
// Compute the mixed model deformation gradient, F_bar
double r1 = pow(fabs(Theta/F.det()),1.0/3.0);
unsigned short int i, j, k;
for (i = 1; i <= nsd; i++)
for (j = 1; j <= nsd; j++)
Fbar(i,j) = r1*F(i,j);
if (nsd == 2) // In 2D we always assume plane strain so set F(3,3)=1
Fbar(3,3) = r1;
else if (nsd != 3)
return false;
// Invert the deformation gradient ==> Fi
double J = F.inverse();
if (J == 0.0) return false;
// Push-forward the basis function gradients to current configuration
dNdx.multiply(dNdX1,F); // dNdx = dNdX * F^-1
#if INT_DEBUG > 0
std::cout <<"NonlinearElasticityULMixed::dNdX ="<< dNdX1;
std::cout <<"NonlinearElasticityULMixed::dNdx ="<< dNdx;
std::cout <<"NonlinearElasticityULMixed::Fbar ="<< Fbar;
#endif
if (eM)
// Integrate the mass matrix
this->formMassMatrix(*eM,N1,X,J*detJW);
if (eS)
// Integrate the load vector due to gravitation and other body forces
this->formBodyForce(*eS,N1,X,J*detJW);
// Evaluate the constitutive relation
F = Fbar;
SymmTensor Eps(3), Sig(3), Sigma(nsd);
double U, Bpres = 0.0, Mpres = 0.0;
if (!this->constitutive(Cmat,Sig,U,Eps=E,X,lHaveStrains))
return false;
#ifdef USE_FTNMAT
// Project the constitutive matrix for the mixed model
int ipsw = INT_DEBUG > 1 ? 9 : 0;
pcst3d_(Cmat.ptr(),Dmat.ptr(),ipsw,6);
#if INT_DEBUG == 1
std::cout <<"NonlinearElasticityULMixed::Dmat ="<< Dmat;
#endif
#else
return false;
#endif
if (lHaveStrains)
{
// Compute mixed strees
Bpres = Sig.trace()/3.0;
Mpres = Press * J/Theta;
Sigma = Sig + (Mpres-Bpres);
// Integrate the geometric stiffness matrix
this->formKG(*eKg,dNdx,Sigma,dVol);
}
#ifdef USE_FTNMAT
mdma3d_(Bpres,Mpres,Sig.ptr(),Dmat.ptr(),ipsw,6);
// Integrate the material stiffness matrix
Dmat *= dVol;
if (nsd == 2)
acckm2d_(N1.size(),dNdx.ptr(),Dmat.ptr(),eKm->ptr());
else
acckm3d_(N1.size(),dNdx.ptr(),Dmat.ptr(),eKm->ptr());
#endif
// Integrate the volumetric change and pressure tangent terms
size_t a, b;
for (a = 1; a <= N1.size(); a++)
for (b = 1; b <= N2.size(); b++)
for (i = 1; i <= nsd; i++)
{
myMats->A[Kut](nsd*(a-1)+i,b) += dNdx(a,i)*N2(b) * Dmat(i,7);
if (nsd == 2)
myMats->A[Kut](nsd*(a-1)+i,b) += dNdx(a,3-i)*N2(b) * Dmat(4,7);
else if (i < 3)
{
j = i + 1;
k = j%3 + 1;
myMats->A[Kut](nsd*(a-1)+i,b) += dNdx(a,3-i)*N2(b) * Dmat(4,7);
myMats->A[Kut](nsd*(a-1)+j,b) += dNdx(a,5-j)*N2(b) * Dmat(5,7);
myMats->A[Kut](nsd*(a-1)+k,b) += dNdx(a,4-k)*N2(b) * Dmat(6,7);
}
myMats->A[Kup](nsd*(a-1)+i,b) += dNdx(a,i)*N2(b) * dVol;
}
myMats->A[Ktt].outer_product(N2,N2*Dmat(7,7),true); // Ktt += N2*N2^T * D77
myMats->A[Ktp].outer_product(N2,N2*(-detJW),true); // Ktp -= N2*N2^T * detJW
if (lHaveStrains)
{
// Compute the small-deformation strain-displacement matrix B from dNdx
if (!this->Elasticity::formBmatrix(dNdx)) return false;
// Integrate the internal forces
Sigma *= -dVol;
if (!Bmat.multiply(Sigma,myMats->b[Ru],true,true)) // Ru -= B^T * sigma
return false;
// Integrate the volumetric change and pressure forces
myMats->b[Rt].add(N2,(Press - Bpres)*detJW); // Rt += N2 * (p-pBar)
myMats->b[Rp].add(N2,(Theta - J)*detJW); // Rp += N2 * (Theta-J)
}
return this->getIntegralResult(elmInt);
}

View File

@@ -0,0 +1,108 @@
// $Id: NonlinearElasticityULMixed.h,v 1.2 2010-12-30 15:25:14 kmo Exp $
//==============================================================================
//!
//! \file NonlinearElasticityULMixed.h
//!
//! \date Dec 22 2010
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Integrand implementations for mixed nonlinear elasticity problems.
//!
//==============================================================================
#ifndef _NONLINEAR_ELASTICITY_UL_MIXED_H
#define _NONLINEAR_ELASTICITY_UL_MIXED_H
#include "NonlinearElasticityUL.h"
#include "ElmMats.h"
/*!
\brief Class representing the integrand of the nonlinear elasticity problem.
\details This class implements a mixed Updated Lagrangian formulation with
continuous pressure and volumetric change fields.
*/
class NonlinearElasticityULMixed : public NonlinearElasticityUL
{
//! \brief Class representing the element matrices of the mixed formulation.
class MixedElmMats : public ElmMats
{
public:
//! \brief Default constructor.
MixedElmMats();
//! \brief Empty destructor.
virtual ~MixedElmMats() {}
//! \brief Returns the element-level Newton matrix.
virtual const Matrix& getNewtonMatrix() const;
//! \brief Returns the element-level right-hand-side vector
//! associated with the Newton matrix.
virtual const Vector& getRHSVector() const;
};
public:
//! \brief The default constructor invokes the parent class constructor.
//! \param[in] n Number of spatial dimensions
//! \param[in] mver Material version parameter
NonlinearElasticityULMixed(unsigned short int n = 3, int mver = 0);
//! \brief Empty destructor.
virtual ~NonlinearElasticityULMixed() {}
//! \brief Prints out problem definition to the given output stream.
virtual void print(std::ostream& os) const;
//! \brief Defines the solution mode before the element assembly is started.
//! \param[in] mode The solution mode to use
virtual void setMode(SIM::SolutionMode mode);
//! \brief Initializes current element for numerical integration.
//! \param[in] MNPC1 Nodal point correspondance for the basis 1
//! \param[in] MNPC2 Nodal point correspondance for the basis 2
//! \param[in] n1 Number of nodes in basis 1 on this patch
virtual bool initElement(const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2, size_t n1);
//! \brief Initializes current element for numerical boundary integration.
//! \param[in] MNPC1 Nodal point correspondance for the basis 1
//! \param[in] MNPC2 Nodal point correspondance for the basis 2
virtual bool initElementBou(const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2, size_t);
//! \brief Evaluates the mixed field problem integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N1 Basis function values, field 1 (displacements)
//! \param[in] N2 Basis function values, field 2 (pressure and volume change)
//! \param[in] dN1dX Basis function gradients, field 1
//! \param[in] dN2dX Basis function gradients, field 2
//! \param[in] X Cartesian coordinates of current integration point
virtual bool evalInt(LocalIntegral*& elmInt, double detJW,
const Vector& N1, const Vector& N2,
const Matrix& dN1dX, const Matrix& dN2dX,
const Vec3& X) const;
//! \brief Evaluates the integrand at a boundary point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N1 Basis function values, field 1
//! \param[in] dN1dX Basis function gradients, field 1
//! \param[in] X Cartesian coordinates of current integration point
//! \param[in] normal Boundary normal vector at current integration point
//!
//! \details The boundary integral is the same as that of the parent class.
//! It does not depend on the pressure and volumetric change fields.
//! Thus, this call is forwarded to the single-field parent class method.
virtual bool evalBou(LocalIntegral*& elmInt, double detJW,
const Vector& N1, const Vector&,
const Matrix& dN1dX, const Matrix&,
const Vec3& X, const Vec3& normal) const
{
return this->NonlinearElasticityUL::evalBou(elmInt,detJW,N1,dN1dX,X,normal);
}
protected:
mutable Matrix Fbar; //!< Mixed model deformation gradient
mutable Matrix Dmat; //!< Projected mixed constitutive matrix
};
#endif

View File

@@ -0,0 +1,89 @@
// $Id: SIMFiniteDefEl.C,v 1.3 2010-12-29 18:51:53 kmo Exp $
//==============================================================================
//!
//! \file SIMFiniteDefEl.C
//!
//! \date Dec 18 2010
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Solution drivers for NURBS-based finite deformation analysis.
//!
//==============================================================================
#include "SIMFiniteDefEl.h"
#include "NonlinearElasticityULMixed.h"
#include "NonlinearElasticityULMX.h"
#include "NeoHookeElasticity.h"
SIMFiniteDefEl2D::SIMFiniteDefEl2D (int form, bool planeStress,
const std::vector<int>& options)
: SIMLinEl2D(SIM::NONE)
{
int mVER = options.size() > 0 ? options[0] : -1; // material version
int pOrd = options.size() > 1 ? options[1] : 0; // pressure field order
switch (form)
{
case SIM::MIXED_QnQn1:
nf[1] = 2; // continuous volumetric change and pressure fields
myProblem = new NonlinearElasticityULMixed(2,mVER);
break;
case SIM::MIXED_QnPn1:
// Local discontinuous volumetric change and pressure fields
myProblem = new NonlinearElasticityULMX(2,mVER,pOrd);
break;
case SIM::UPDATED_LAGRANGE:
myProblem = new NonlinearElasticityUL(2,planeStress,mVER);
break;
case SIM::TOTAL_LAGRANGE:
myProblem = new NonlinearElasticityTL(2,planeStress);
break;
case SIM::NONLINEAR: // Old tensor-based TL-formulation
myProblem = new NonlinearElasticity(2,planeStress);
break;
default:
std::cerr <<" *** SIMFiniteDefEl2D: Unknown problem formulation "
<< form << std::endl;
}
}
SIMFiniteDefEl3D::SIMFiniteDefEl3D (bool checkRHS, int form,
const std::vector<int>& options)
: SIMLinEl3D(checkRHS,SIM::NONE)
{
int mVER = options.size() > 0 ? options[0] : -1; // material version
int pOrd = options.size() > 1 ? options[1] : 0; // pressure field order
switch (form)
{
case SIM::MIXED_QnQn1:
nf[1] = 2; // continuous volumetric change and pressure fields
myProblem = new NonlinearElasticityULMixed(3,mVER);
break;
case SIM::MIXED_QnPn1:
// Local discontinuous volumetric change and pressure fields
myProblem = new NonlinearElasticityULMX(3,mVER,pOrd);
break;
case SIM::UPDATED_LAGRANGE:
myProblem = new NonlinearElasticityUL(3,false,mVER);
break;
case SIM::TOTAL_LAGRANGE:
myProblem = new NonlinearElasticityTL();
break;
case SIM::NONLINEAR: // Old tensor-based TL-formulation
myProblem = new NonlinearElasticity();
break;
case SIM::NEOHOOKE: // NeoHookean TL-formulation (not working)
myProblem = new NeoHookeElasticity();
break;
case SIM::NEOHOOKE_IV: // NeoHookean isochoric/volumetric TL-formulation
myProblem = new NeoHookeElasticityIV();
break;
default:
std::cerr <<" *** SIMFiniteDefEl3D: Unknown problem formulation "
<< form << std::endl;
}
}

View File

@@ -0,0 +1,71 @@
// $Id: SIMFiniteDefEl.h,v 1.2 2010-12-29 18:51:53 kmo Exp $
//==============================================================================
//!
//! \file SIMFiniteDefEl.h
//!
//! \date Dec 18 2010
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Solution drivers for NURBS-based finite deformation analysis.
//!
//==============================================================================
#ifndef _SIM_FINITE_DEF_EL_H
#define _SIM_FINITE_DEF_EL_H
#include "SIMLinEl2D.h"
#include "SIMLinEl3D.h"
namespace SIM
{
//! \brief Enum defining various finite deformation formulations.
enum NlFormulation
{
TOTAL_LAGRANGE = 3,
NEOHOOKE = 4,
NEOHOOKE_IV = 5,
UPDATED_LAGRANGE = 6,
MIXED_QnPn1 = 7,
MIXED_QnQn1 = 8
};
};
/*!
\brief Driver class for 2D isogeometric finite deformation analysis.
*/
class SIMFiniteDefEl2D : public SIMLinEl2D
{
public:
//! \brief Default constructor.
//! \param[in] form Problem formulation option
//! \param[in] planeStress Plane stress/plane strain option
//! \param[in] options Additional solution formulation options
SIMFiniteDefEl2D(int form = SIM::TOTAL_LAGRANGE, bool planeStress = false,
const std::vector<int>& options = std::vector<int>());
//! \brief Empty destructor.
virtual ~SIMFiniteDefEl2D() {}
};
/*!
\brief Driver class for 3D isogeometric finite deformation analysis.
*/
class SIMFiniteDefEl3D : public SIMLinEl3D
{
public:
//! \brief Default constructor.
//! \param[in] checkRHS If \e true, ensure the model is in a right-hand system
//! \param[in] form Problem formulation option
//! \param[in] options Additional solution formulation options
SIMFiniteDefEl3D(bool checkRHS = false, int form = SIM::TOTAL_LAGRANGE,
const std::vector<int>& options = std::vector<int>());
//! \brief Empty destructor.
virtual ~SIMFiniteDefEl3D() {}
};
#endif

View File

@@ -0,0 +1,135 @@
C $Id: accKM.f,v 1.1 2010-12-29 18:51:53 kmo Exp $
C=======================================================================
C!
C! \file accKM.f
C!
C! \date Dec 18 2010
C!
C! \author Knut Morten Okstad / SINTEF
C!
C! \brief Material stiffness accumulation routines.
C!
C=======================================================================
C
subroutine accKM2D (NENOD,Shpf,D,EKt)
C
implicit none
C
integer NENOD
double precision Shpf(NENOD,2), D(7,7)
double precision EKt(2,NENOD,2,NENOD)
C
integer a, b, i, j
double precision BBC(2,4)
C
do a = 1, NENOD
C
do j = 1, 4
BBC(1,j) = Shpf(a,1)*D(1,j)
& + Shpf(a,2)*D(4,j)
BBC(2,j) = Shpf(a,1)*D(4,j)
& + Shpf(a,2)*D(2,j)
end do
C
do b = 1, NENOD
C
EKt(1,a,1,b) = EKt(1,a,1,b)
& + BBC(1,1) * Shpf(b,1)
& + BBC(1,4) * Shpf(b,2)
C
EKt(1,a,2,b) = EKt(1,a,2,b)
& + BBC(1,4) * Shpf(b,1)
& + BBC(1,2) * Shpf(b,2)
C
EKt(2,a,1,b) = EKt(2,a,1,b)
& + BBC(2,1) * Shpf(b,1)
& + BBC(2,4) * Shpf(b,2)
C
EKt(2,a,2,b) = EKt(2,a,2,b)
& + BBC(2,4) * Shpf(b,1)
& + BBC(2,2) * Shpf(b,2)
C
end do
C
end do
C
end
C
C
subroutine accKM3D (NENOD,Shpf,D,EKt)
C
implicit none
C
integer NENOD
double precision Shpf(NENOD,3), D(7,7)
double precision EKt(3,NENOD,3,NENOD)
C
integer a, b, i, j
double precision BBC(3,6)
C
do a = 1, NENOD
C
do j = 1, 6
BBC(1,j) = Shpf(a,1)*D(1,j)
& + Shpf(a,2)*D(4,j)
& + Shpf(a,3)*D(6,j)
BBC(2,j) = Shpf(a,1)*D(4,j)
& + Shpf(a,2)*D(2,j)
& + Shpf(a,3)*D(5,j)
BBC(3,j) = Shpf(a,1)*D(6,j)
& + Shpf(a,2)*D(5,j)
& + Shpf(a,3)*D(3,j)
end do
C
do b = 1, NENOD
C
EKt(1,a,1,b) = EKt(1,a,1,b)
& + BBC(1,1) * Shpf(b,1)
& + BBC(1,4) * Shpf(b,2)
& + BBC(1,6) * Shpf(b,3)
C
EKt(1,a,2,b) = EKt(1,a,2,b)
& + BBC(1,4) * Shpf(b,1)
& + BBC(1,2) * Shpf(b,2)
& + BBC(1,5) * Shpf(b,3)
C
EKt(1,a,3,b) = EKt(1,a,3,b)
& + BBC(1,6) * Shpf(b,1)
& + BBC(1,5) * Shpf(b,2)
& + BBC(1,3) * Shpf(b,3)
C
EKt(2,a,1,b) = EKt(2,a,1,b)
& + BBC(2,1) * Shpf(b,1)
& + BBC(2,4) * Shpf(b,2)
& + BBC(2,6) * Shpf(b,3)
C
EKt(2,a,2,b) = EKt(2,a,2,b)
& + BBC(2,4) * Shpf(b,1)
& + BBC(2,2) * Shpf(b,2)
& + BBC(2,5) * Shpf(b,3)
C
EKt(2,a,3,b) = EKt(2,a,3,b)
& + BBC(2,6) * Shpf(b,1)
& + BBC(2,5) * Shpf(b,2)
& + BBC(2,3) * Shpf(b,3)
C
EKt(3,a,1,b) = EKt(3,a,1,b)
& + BBC(3,1) * Shpf(b,1)
& + BBC(3,4) * Shpf(b,2)
& + BBC(3,6) * Shpf(b,3)
C
EKt(3,a,2,b) = EKt(3,a,2,b)
& + BBC(3,4) * Shpf(b,1)
& + BBC(3,2) * Shpf(b,2)
& + BBC(3,5) * Shpf(b,3)
C
EKt(3,a,3,b) = EKt(3,a,3,b)
& + BBC(3,6) * Shpf(b,1)
& + BBC(3,5) * Shpf(b,2)
& + BBC(3,3) * Shpf(b,3)
C
end do
C
end do
C
end

View File

@@ -0,0 +1,153 @@
C $Id: accKMx.f,v 1.2 2010-12-20 17:01:21 kmo Exp $
C=======================================================================
C!
C! \file accKMx.f
C!
C! \date Dec 18 2010
C!
C! \author Knut Morten Okstad / SINTEF
C!
C! \brief Material stiffness accumulation routines for B-bar methods.
C!
C=======================================================================
C
subroutine accKMx2D (NENOD,Shpf,Shpbar,D,EKt)
C
implicit none
C
integer NENOD
double precision Shpf(NENOD,2), Shpbar(NENOD,2), D(7,7)
double precision EKt(2,NENOD,2,NENOD)
C
integer a, b, i, j
double precision BBC(2,7)
C
do a = 1, NENOD
C
do j = 1, 7
BBC(1,j) = Shpf(a,1)*D(1,j)
& + Shpf(a,2)*D(4,j)
& + Shpbar(a,1)*D(7,j)
BBC(2,j) = Shpf(a,1)*D(4,j)
& + Shpf(a,2)*D(2,j)
& + Shpbar(a,2)*D(7,j)
end do
C
do b = 1, NENOD
C
EKt(1,a,1,b) = EKt(1,a,1,b)
& + BBC(1,1) * Shpf(b,1)
& + BBC(1,4) * Shpf(b,2)
& + BBC(1,7) * Shpbar(b,1)
C
EKt(1,a,2,b) = EKt(1,a,2,b)
& + BBC(1,4) * Shpf(b,1)
& + BBC(1,2) * Shpf(b,2)
& + BBC(1,7) * Shpbar(b,2)
C
EKt(2,a,1,b) = EKt(2,a,1,b)
& + BBC(2,1) * Shpf(b,1)
& + BBC(2,4) * Shpf(b,2)
& + BBC(2,7) * Shpbar(b,1)
C
EKt(2,a,2,b) = EKt(2,a,2,b)
& + BBC(2,4) * Shpf(b,1)
& + BBC(2,2) * Shpf(b,2)
& + BBC(2,7) * Shpbar(b,2)
C
end do
C
end do
C
end
C
C
subroutine accKMx3D (NENOD,Shpf,Shpbar,D,EKt)
C
implicit none
C
integer NENOD
double precision Shpf(NENOD,3), Shpbar(NENOD,3), D(7,7)
double precision EKt(3,NENOD,3,NENOD)
C
integer a, b, i, j
double precision BBC(3,7)
C
do a = 1, NENOD
C
do j = 1, 7
BBC(1,j) = Shpf(a,1)*D(1,j)
& + Shpf(a,2)*D(4,j)
& + Shpf(a,3)*D(6,j)
& + Shpbar(a,1)*D(7,j)
BBC(2,j) = Shpf(a,1)*D(4,j)
& + Shpf(a,2)*D(2,j)
& + Shpf(a,3)*D(5,j)
& + Shpbar(a,2)*D(7,j)
BBC(3,j) = Shpf(a,1)*D(6,j)
& + Shpf(a,2)*D(5,j)
& + Shpf(a,3)*D(3,j)
& + Shpbar(a,3)*D(7,j)
end do
C
do b = 1, NENOD
C
EKt(1,a,1,b) = EKt(1,a,1,b)
& + BBC(1,1) * Shpf(b,1)
& + BBC(1,4) * Shpf(b,2)
& + BBC(1,6) * Shpf(b,3)
& + BBC(1,7) * Shpbar(b,1)
C
EKt(1,a,2,b) = EKt(1,a,2,b)
& + BBC(1,4) * Shpf(b,1)
& + BBC(1,2) * Shpf(b,2)
& + BBC(1,5) * Shpf(b,3)
& + BBC(1,7) * Shpbar(b,2)
C
EKt(1,a,3,b) = EKt(1,a,3,b)
& + BBC(1,6) * Shpf(b,1)
& + BBC(1,5) * Shpf(b,2)
& + BBC(1,3) * Shpf(b,3)
& + BBC(1,7) * Shpbar(b,3)
C
EKt(2,a,1,b) = EKt(2,a,1,b)
& + BBC(2,1) * Shpf(b,1)
& + BBC(2,4) * Shpf(b,2)
& + BBC(2,6) * Shpf(b,3)
& + BBC(2,7) * Shpbar(b,1)
C
EKt(2,a,2,b) = EKt(2,a,2,b)
& + BBC(2,4) * Shpf(b,1)
& + BBC(2,2) * Shpf(b,2)
& + BBC(2,5) * Shpf(b,3)
& + BBC(2,7) * Shpbar(b,2)
C
EKt(2,a,3,b) = EKt(2,a,3,b)
& + BBC(2,6) * Shpf(b,1)
& + BBC(2,5) * Shpf(b,2)
& + BBC(2,3) * Shpf(b,3)
& + BBC(2,7) * Shpbar(b,3)
C
EKt(3,a,1,b) = EKt(3,a,1,b)
& + BBC(3,1) * Shpf(b,1)
& + BBC(3,4) * Shpf(b,2)
& + BBC(3,6) * Shpf(b,3)
& + BBC(3,7) * Shpbar(b,1)
C
EKt(3,a,2,b) = EKt(3,a,2,b)
& + BBC(3,4) * Shpf(b,1)
& + BBC(3,2) * Shpf(b,2)
& + BBC(3,5) * Shpf(b,3)
& + BBC(3,7) * Shpbar(b,2)
C
EKt(3,a,3,b) = EKt(3,a,3,b)
& + BBC(3,6) * Shpf(b,1)
& + BBC(3,5) * Shpf(b,2)
& + BBC(3,3) * Shpf(b,3)
& + BBC(3,7) * Shpbar(b,3)
C
end do
C
end do
C
end

View File

@@ -0,0 +1,313 @@
// $Id: main_NonLinEl.C,v 1.6 2011-02-08 09:32:18 kmo Exp $
//==============================================================================
//!
//! \file main_NonLinEl.C
//!
//! \date Jun 1 2010
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Main program for the isogeometric nonlinear elasticity solver.
//!
//==============================================================================
#include "NonLinSIM.h"
#include "SIMFiniteDefEl.h"
#include "LinAlgInit.h"
#include "Profiler.h"
#include "VTF.h"
#include <fstream>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
/*!
\brief Main program for the isogeometric nonlinear elasticity solver.
The input to the program is specified through the following
command-line arguments. The arguments may be given in arbitrary order.
\arg \a input-file : Input file with model definition
\arg -dense : Use the dense LAPACK matrix equation solver
\arg -spr : Use the SPR direct equation solver
\arg -superlu : Use the sparse SuperLU equation solver
\arg -samg : Use the sparse algebraic multi-grid equation solver
\arg -petsc : Use equation solver from PETSc library
\arg -nGauss \a n : Number of Gauss points over a knot-span in each direction
\arg -vtf \a format : VTF-file format (0=ASCII, 1=BINARY)
\arg -nviz \a nviz : Number of visualization points over each knot-span
\arg -nu \a nu : Number of visualization points per knot-span in u-direction
\arg -nv \a nv : Number of visualization points per knot-span in v-direction
\arg -nw \a nw : Number of visualization points per knot-span in w-direction
\arg -skip2nd : Skip VTF-output of secondary solution fields
\arg -saveInc \a dtSave : Time increment between each result save to VTF
\arg -dumpInc \a dtDump [raw] : Time increment between each solution dump
\arg -ignore \a p1, \a p2, ... : Ignore these patches in the analysis
\arg -check : Data check only, read model and output to VTF (no solution)
\arg -checkRHS : Check that the patches are modelled in a right-hand system
\arg -fixDup : Resolve co-located nodes by merging them into a single node
\arg -2D : Use two-parametric simulation driver
\arg -Tensor : Use tensorial total Lagrangian formulation (slow...)
\arg -UL \a mver : Use updated Lagrangian formulation with nonlinear material
\arg -MX \a pord : Mixed formulation with internal discontinuous pressure
\arg -mixed : Mixed formulation with continuous pressure and volumetric change
\arg -lag : Use Lagrangian basis functions instead of splines/NURBS
\arg -spec : Use Spectral basis functions instead of splines/NURBS
*/
int main (int argc, char** argv)
{
Profiler prof(argv[0]);
SystemMatrix::Type solver = SystemMatrix::SPARSE;
int form = SIM::TOTAL_LAGRANGE;
int nGauss = 4;
int format = 0;
int n[3] = { 2, 2, 2 };
std::vector<int> ignoredPatches;
double dtSave = 0.0;
double dtDump = 0.0;
bool dumpWithID = true;
bool skip2nd = false;
bool checkRHS = false;
bool fixDup = false;
char twoD = false;
char* infile = 0;
LinAlgInit linalg(argc,argv);
std::vector<int> options(2,-1);
for (int i = 1; i < argc; i++)
if (!strcmp(argv[i],"-dense"))
solver = SystemMatrix::DENSE;
else if (!strcmp(argv[i],"-spr"))
solver = SystemMatrix::SPR;
else if (!strncmp(argv[i],"-superlu",8))
{
solver = SystemMatrix::SPARSE;
if (isdigit(argv[i][8]))
SIMbase::num_threads_SLU = atoi(argv[i]+8);
}
else if (!strcmp(argv[i],"-samg"))
solver = SystemMatrix::SAMG;
else if (!strcmp(argv[i],"-petsc"))
solver = SystemMatrix::PETSC;
else if (!strcmp(argv[i],"-nGauss") && i < argc-1)
nGauss = atoi(argv[++i]);
else if (!strcmp(argv[i],"-vtf") && i < argc-1)
format = atoi(argv[++i]);
else if (!strcmp(argv[i],"-nviz") && i < argc-1)
n[0] = n[1] = n[2] = atoi(argv[++i]);
else if (!strcmp(argv[i],"-nu") && i < argc-1)
n[0] = atoi(argv[++i]);
else if (!strcmp(argv[i],"-nv") && i < argc-1)
n[1] = atoi(argv[++i]);
else if (!strcmp(argv[i],"-nw") && i < argc-1)
n[2] = atoi(argv[++i]);
else if (!strcmp(argv[i],"-skip2nd"))
skip2nd = true;
else if (!strcmp(argv[i],"-saveInc") && i < argc-1)
dtSave = atof(argv[++i]);
else if (!strcmp(argv[i],"-dumpInc") && i < argc-1)
{
dtDump = atof(argv[++i]);
if (++i < argc && !strcmp(argv[i],"raw"))
dumpWithID = false;
else
--i;
}
else if (!strcmp(argv[i],"-ignore"))
while (i < argc-1 && isdigit(argv[i+1][0]))
{
char* endp = 0;
int endVal = 0;
ignoredPatches.push_back(strtol(argv[++i],&endp,10));
if (endp && *endp == ':')
endVal = strtol(endp+1,&endp,10);
while (ignoredPatches.back() < endVal)
ignoredPatches.push_back(ignoredPatches.back()+1);
}
else if (!strcmp(argv[i],"-vox") && i < argc-1)
VTF::vecOffset[0] = atof(argv[++i]);
else if (!strcmp(argv[i],"-voy") && i < argc-1)
VTF::vecOffset[1] = atof(argv[++i]);
else if (!strcmp(argv[i],"-voz") && i < argc-1)
VTF::vecOffset[2] = atof(argv[++i]);
else if (!strcmp(argv[i],"-check"))
form += 100;
else if (!strcmp(argv[i],"-checkRHS"))
checkRHS = true;
else if (!strcmp(argv[i],"-fixDup"))
fixDup = true;
else if (!strncmp(argv[i],"-2D",3))
twoD = strcmp(argv[i],"-2Dpstrain") ? 1 : 2;
else if (!strcmp(argv[i],"-UL"))
{
if (form < SIM::UPDATED_LAGRANGE)
form = SIM::UPDATED_LAGRANGE;
if (i < argc-1 && isdigit(argv[i+1][0]))
options[0] = atoi(argv[++i]);
}
else if (!strcmp(argv[i],"-MX"))
{
form = SIM::MIXED_QnPn1;
if (i < argc-1 && isdigit(argv[i+1][0]))
options[1] = atoi(argv[++i]);
else
options[1] = 0;
}
else if (!strcmp(argv[i],"-mixed"))
form = SIM::MIXED_QnQn1;
else if (!strcmp(argv[i],"-Tensor"))
form = SIM::NONLINEAR;
else if (!strcmp(argv[i],"-NeoHooke"))
form = SIM::NEOHOOKE;
else if (!strcmp(argv[i],"-NeoHookeIV"))
form = SIM::NEOHOOKE_IV;
else if (!strncmp(argv[i],"-lag",4))
SIMbase::discretization = SIMbase::Lagrange;
else if (!strncmp(argv[i],"-spec",5))
SIMbase::discretization = SIMbase::Spectral;
else if (!infile)
infile = argv[i];
else
std::cerr <<" ** Unknown option ignored: "<< argv[i] << std::endl;
if (!infile)
{
std::cout <<"usage: "<< argv[0]
<<" <inputfile> [-dense|-spr|-superlu[<nt>]|-samg|-petsc]\n "
<<" [-lag] [-spec] [-2D[pstrain]] [-UL [<mVER>]] [-MX [<p>]|-mixe"
<<"d] [-nGauss <n>]\n [-vtf <format>] [-nviz <nviz>]"
<<" [-nu <nu>] [-nv <nv>] [-nw <nw>]\n "
<<" [-saveInc <dtSave>] [-dumpInc <dtDump> [raw]]\n "
<<" [-ignore <p1> <p2> ...] [-fixDup] [-checkRHS] [-check]\n";
return 0;
}
if (linalg.myPid == 0)
std::cout <<"\n >>> Spline FEM Nonlinear Elasticity solver <<<"
<<"\n ==============================================\n"
<<"\nInput file: "<< infile
<<"\nEquation solver: "<< solver
<<"\nNumber of Gauss points: "<< nGauss
<<"\nVTF file format: "<< (format ? "BINARY":"ASCII")
<<"\nNumber of visualization points: "
<< n[0] <<" "<< n[1];
if (twoD)
n[2] = 1;
else if (linalg.myPid == 0)
std::cout <<" "<< n[2];
if (linalg.myPid == 0)
{
if (dtSave > 0.0)
std::cout <<"\nTime between each result save: "<< dtSave;
if (dtDump > 0.0)
std::cout <<"\nTime between each primary solution dump: "<< dtDump;
if (SIMbase::discretization == SIMbase::Lagrange)
std::cout <<"\nLagrangian basis functions are used";
else if (SIMbase::discretization == SIMbase::Spectral)
std::cout <<"\nSpectral basis functions are used";
if (fixDup)
std::cout <<"\nCo-located nodes will be merged";
if (checkRHS)
std::cout <<"\nCheck that each patch has a right-hand coordinate system";
if (!ignoredPatches.empty())
{
std::cout <<"\nIgnored patches:";
for (size_t i = 0; i < ignoredPatches.size(); i++)
std::cout <<" "<< ignoredPatches[i];
}
std::cout << std::endl;
}
utl::profiler->start("Model input");
// Create the finite deformation elasticity model
SIMbase* model;
if (twoD)
model = new SIMFiniteDefEl2D(form%100,twoD==1,options);
else
model = new SIMFiniteDefEl3D(checkRHS,form%100,options);
// Read in solver and model definitions
NonLinSIM simulator(model);
if (!simulator.read(infile))
return 1;
utl::profiler->stop("Model input");
// Preprocess the model and establish data structures for the algebraic system
model->printProblem(std::cout);
if (!model->preprocess(ignoredPatches,fixDup))
return 2;
// Save FE model to VTF file for visualization
if (!simulator.saveModel(infile,format,n))
return 3;
std::ostream* oss = 0;
if (dtDump > 0.0 && !dumpWithID)
{
// Write (refined?) model to g2-file
strcat(strtok(infile,"."),".g2");
std::cout <<"\nWriting updated g2-file "<< infile << std::endl;
std::ofstream osg(infile);
model->dumpGeometry(osg);
// Open ASCII file for solution dump
strcat(strtok(infile,"."),".sol");
oss = new std::ofstream(infile);
*oss <<"#NPoints="<< model->getNoNodes() <<"\n";
}
// Define the initial configuration
NonLinSIM::SolvePrm params;
simulator.init(params);
const double epsT = 1.0e-6;
if (dtDump <= 0.0) dtDump = params.stopTime + 1.0;
double nextDump = params.time.t + dtDump;
double nextSave = params.time.t + dtSave;
int iStep = 0; // Save initial state to VTF
if (params.multiSteps())
if (!simulator.saveStep(-(++iStep),params.time.t,n,skip2nd))
return 4;
if (form >= 100)
return 0; // model check
// Initialize the linear solver
model->initSystem(solver,1,1);
model->setAssociatedRHS(0,0);
model->setQuadratureRule(nGauss);
// Invoke the time/load-step loop
while (simulator.advanceStep(params))
{
// Solve the nonlinear FE problem at this load step
if (!simulator.solveStep(params,SIM::STATIC))
return 5;
// Dump primary solution for inspection or external processing
if (params.time.t + epsT*params.time.dt > nextDump)
{
if (dumpWithID)
simulator.dumpStep(params.step,params.time.t,std::cout);
else
simulator.dumpStep(params.step,params.time.t,*oss,false);
nextDump = params.time.t + dtDump;
}
// Save solution variables to VTF for visualization
if (params.time.t + epsT*params.time.dt > nextSave)
if (simulator.saveStep(++iStep,params.time.t,n,skip2nd))
nextSave = params.time.t + dtSave;
else
return 6;
}
if (oss) delete oss;
return 0;
}

View File

@@ -0,0 +1,205 @@
C $Id: stiff_TL.f,v 1.1 2010-12-18 20:03:54 kmo Exp $
C=======================================================================
C!
C! \file stiff_TL.f
C!
C! \date Jun 8 2010
C!
C! \author Knut Morten Okstad / SINTEF
C!
C! \brief Routines for Total Lagrangian material stiffness integration.
C!
C=======================================================================
C
subroutine stiff_TL2D (NENOD,detJW,dNdX,F,Cmat,EK)
C
implicit none
C
integer NENOD, a, b, i, j, k, l, m, n
double precision detJW, dNdX(NENOD,2), F(2,2), Cmat(3,3)
double precision EK(2,NENOD,2,NENOD), C(2,2,2,2), km
C
do 200 i = 1, 2
do 200 j = i, 2
if (i == j) then
a = i
else
a = 3
end if
do 100 k = 1, 2
do 100 l = k, 2
if (k == l) then
b = k
else
b = 3
end if
C(i,j,k,l) = Cmat(a,b)
C(j,i,k,l) = Cmat(a,b)
C(i,j,l,k) = Cmat(a,b)
C(j,i,l,k) = Cmat(a,b)
100 continue
200 continue
C
do 450 a = 1, NENOD
do 450 b = a, NENOD
do 400 m = 1, 2
do 400 n = 1, 2
km = 0.0D0
do 350 i = 1, 2
do 350 j = 1, 2
do 300 k = 1, 2
do 300 l = 1, 2
km = km + dNdX(a,i)*F(m,j)*C(i,j,k,l)*F(n,k)*dNdX(b,l)
300 continue
350 continue
EK(m,a,n,b) = EK(m,a,n,b) + km*detJW
if (b > a) then
EK(n,b,m,a) = EK(n,b,m,a) + km*detJW
end if
400 continue
450 continue
C
end
C
C
subroutine stiff_TL3D (NENOD,detJW,dNdX,F,Cmat,EK)
C
implicit none
C
integer NENOD, a, b, i, j, k, l, m, n
double precision detJW, DNDX(NENOD,3), F(3,3), Cmat(6,6)
double precision EK(3,NENOD,3,NENOD), C(3,3,3,3), km
C
do 200 i = 1, 3
do 200 j = i, 3
if (i == j) then
a = i
else if (i == 1) then
a = j+2
else
a = 6
end if
do 100 k = 1, 3
do 100 l = k, 3
if (k == l) then
b = k
else if (k == 1) then
b = l+2
else
b = 6
end if
C(i,j,k,l) = Cmat(a,b)
C(j,i,k,l) = Cmat(a,b)
C(i,j,l,k) = Cmat(a,b)
C(j,i,l,k) = Cmat(a,b)
100 continue
200 continue
C
do 450 a = 1, NENOD
do 450 b = a, NENOD
do 400 m = 1, 3
do 400 n = 1, 3
km = 0.0D0
do 350 i = 1, 3
do 350 j = 1, 3
do 300 k = 1, 3
do 300 l = 1, 3
km = km + dNdX(a,i)*F(m,j)*C(i,j,k,l)*F(n,k)*dNdX(b,l)
300 continue
350 continue
EK(m,a,n,b) = EK(m,a,n,b) + km*detJW
if (b > a) then
EK(n,b,m,a) = EK(n,b,m,a) + km*detJW
end if
400 continue
450 continue
C
end
C
C
subroutine stiff_TL2D_isoel (NENOD,detJW,dNdX,F,C1,C2,C3,EK)
C
C Use this subroutine in 2D, but only for isotrophic materials,
C where the fourth-order constitutive tensor is assumed to only
C contain three different values, C1, C2, C3 (and zeroes):
C
C C1 = C(1,1,1,1) = C(2,2,2,2)
C C2 = C(1,1,2,2) = C(2,2,1,1
C C3 = C(1,2,1,2) = C(1,2,2,1) = C(2,1,1,2) = C(2,1,2,1)
C
implicit none
C
integer NENOD, a, b, m, n
double precision detJW, DNDX(NENOD,2), F(2,2), C1, C2, C3
double precision EK(2,NENOD,2,NENOD), ek1, ek2, ek3
C
do 350 a = 1, NENOD
do 350 b = a, NENOD
do 300 m = 1, 2
do 300 n = 1, 2
ek1 = 0.0D0
ek2 = 0.0D0
ek3 = 0.0D0
ek1 = dNdX(a,1)*F(m,1)*F(n,1)*dNdX(b,1)
+ + dNdX(a,2)*F(m,2)*F(n,2)*dNdX(b,2)
ek2 = dNdX(a,1)*F(m,1)*F(n,2)*dNdX(b,2)
+ + dNdX(a,2)*F(m,2)*F(n,1)*dNdX(b,1)
ek3 = dNdX(a,1)*F(m,2)*(F(n,1)*dNdX(b,2) +
+ F(n,2)*dNdX(b,1))
+ + dNdX(a,2)*F(m,1)*(F(n,2)*dNdX(b,1) +
+ F(n,1)*dNdX(b,2))
EK(m,a,n,b) = EK(m,a,n,b) + (ek1*C1+ek2*C2+ek3*C3)*detJW
if (b > a) then
EK(n,b,m,a) = EK(n,b,m,a) + (ek1*C1+ek2*C2+ek3*C3)*detJW
end if
300 continue
350 continue
C
end
C
C
subroutine stiff_TL3D_isoel (NENOD,detJW,dNdX,F,C1,C2,C3,EK)
C
C Use this subroutine in 3D, but only for isotrophic materials,
C where the fourth-order constitutive tensor is assumed to only
C contain three different values, C1, C2, C3 (and zeroes):
C
C C1 = C(1,1,1,1) = C(2,2,2,2) = C(3,3,3,3)
C C2 = C(1,1,2,2) = C(1,1,3,3) = C(2,2,3,3)
C = C(2,2,1,1) = C(3,3,1,1) = C(3,3,2,2)
C C3 = C(1,2,1,2) = C(1,2,2,1) = C(2,1,1,2) = C(2,1,2,1)
C = C(1,3,1,3) = C(1,3,3,1) = C(3,1,1,3) = C(3,1,3,1)
C = C(2,3,2,3) = C(2,3,3,2) = C(3,2,2,3) = C(3,2,3,2)
C
implicit none
C
integer NENOD, a, b, i, j, m, n
double precision detJW, DNDX(NENOD,3), F(3,3), C1, C2, C3
double precision EK(3,NENOD,3,NENOD), ek1, ek2, ek3
C
do 350 a = 1, NENOD
do 350 b = a, NENOD
do 300 m = 1, 3
do 300 n = 1, 3
ek1 = 0.0D0
ek2 = 0.0D0
ek3 = 0.0D0
do 200 i = 1, 3
ek1 = ek1 + dNdX(a,i)*F(m,i)*F(n,i)*dNdX(b,i)
do 100 j = i+1, 3
ek2 = ek2 + dNdX(a,i)*F(m,i)*F(n,j)*dNdX(b,j)
+ + dNdX(a,j)*F(m,j)*F(n,i)*dNdX(b,i)
ek3 = ek3 + dNdX(a,i)*F(m,j)*(F(n,i)*dNdX(b,j) +
+ F(n,j)*dNdX(b,i))
+ + dNdX(a,j)*F(m,i)*(F(n,j)*dNdX(b,i) +
+ F(n,i)*dNdX(b,j))
100 continue
200 continue
EK(m,a,n,b) = EK(m,a,n,b) + (ek1*C1+ek2*C2+ek3*C3)*detJW
if (b > a) then
EK(n,b,m,a) = EK(n,b,m,a) + (ek1*C1+ek2*C2+ek3*C3)*detJW
end if
300 continue
350 continue
C
end

76
Apps/Makefile Normal file
View File

@@ -0,0 +1,76 @@
# $Id: Makefile,v 1.13 2011-02-04 17:31:40 kmo Exp $
#===============================================================================
#
# File: Makefile
#
# Created: Jan 12 2010
#
# Author: Knut Morten Okstad / SINTEF
#
# Description:
# Makefile for spline-based isogeometric finite element simulation applications.
#
#===============================================================================
# Comment in this line if PETSc and SLEPc is used
#include ${SLEPC_DIR}/conf/slepc_common_variables
SRCS0 = SIMLinEl2D.C SIMLinEl3D.C \
FiniteDefElasticity/NonlinearElasticity.C \
FiniteDefElasticity/NonlinearElasticityTL.C
SRCS1 = main_LinEl3D.C $(SRCS0)
SRCS2 = main_LinElStep.C $(SRCS0)
OBJS1 = $(SRCS1:.C=.o) FiniteDefElasticity/stiff_TL.o
OBJS2 = $(SRCS2:.C=.o) FiniteDefElasticity/stiff_TL.o
APP1 = linel3D
APP2 = linelStep
MYLIB = ../libSplineFEM.a
GOTLIB = -L../GoTools/lib -lGoTrivariate -lGoToolsCore
BLASLIB = -larpack -llapack -lblas
SLULIB = -lsuperlu
#SLULIB = -L../SuperLU_MT/lib -lsuperlu -lpthread
#SAMLIB = -L../SAM -lSPR -lASM -lEIG -lMAT -lMDR -lUTL
#SAMGLIB = -L../SAMG/samg -lamg_coo
# Comment in this line if PETSc and SLEPc is used
#PETSCOPT = -DHAS_PETSC $(CC_INCLUDES)
# Comment in this line if parallel PETSc and SLEPc are used
#PETSCOPT = -DHAS_PETSC -DPARALLEL_PETSC $(CC_INCLUDES)
VTFLIB = -L../GLviewExpressWriter/lib -lVTFExpressAPI
LIBS = $(MYLIB) $(GOTLIB) $(VTFLIB) \
$(SAMLIB) $(SAMGLIB) $(SLULIB) $(SLEPC_LIB) $(BLASLIB)
FC = gfortran
CXX = g++
#CDEBUG = -g -DINDEX_CHECK
CDEBUG = -O2
FFLAGS = -O3 -funroll-loops
CXXFLAGS = $(CDEBUG) -Dreal=double -DUSE_CBLAS \
-I. -I../Integrands -I../SIM -I../ASM -I../LinAlg -I../Utility \
-IFiniteDefElasticity $(PETSCOPT)
all: $(APP1) $(APP2)
$(APP1): $(OBJS1) $(MYLIB) Makefile
$(CXX) -o $@ $(OBJS1) $(LIBS) -lgfortran
$(APP2): $(OBJS2) $(MYLIB) Makefile
$(CXX) -o $@ $(OBJS2) $(LIBS) -lgfortran
../libSplineFEM.a:
cd ..; make
clean:
rm -f $(OBJS1) $(OBJS2) $(OBJS3)
sed '1,/^# DO NOT DELETE/w Makefile.tmp' < Makefile > /dev/null
mv Makefile.tmp Makefile
depend:
sed '1,/^# DO NOT DELETE/w Makefile.tmp' < Makefile > /dev/null
$(CXX) -M $(CXXFLAGS) *.C >> Makefile.tmp
mv Makefile Makefile.bak
mv Makefile.tmp Makefile
# DO NOT DELETE

293
Apps/SIMLinEl2D.C Normal file
View File

@@ -0,0 +1,293 @@
// $Id: SIMLinEl2D.C,v 1.18 2011-02-08 09:06:02 kmo Exp $
//==============================================================================
//!
//! \file SIMLinEl2D.C
//!
//! \date Dec 04 2010
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Solution driver for 2D NURBS-based linear elastic FEM analysis.
//!
//==============================================================================
#include "SIMLinEl2D.h"
#include "LinearElasticity.h"
#include "NonlinearElasticity.h"
#include "AnalyticSolutions.h"
#include "Functions.h"
#include "Utilities.h"
#include "Vec3Oper.h"
#include "Property.h"
#include "Tensor.h"
#include <string.h>
SIMLinEl2D::SIMLinEl2D (int form, bool planeStress)
{
switch (form)
{
case SIM::LINEAR:
myProblem = new LinearElasticity(2,planeStress);
break;
case SIM::NONLINEAR:
myProblem = new NonlinearElasticity(2,planeStress);
break;
}
asol = 0;
}
SIMLinEl2D::~SIMLinEl2D ()
{
if (asol) delete asol;
}
bool SIMLinEl2D::parse (char* keyWord, std::istream& is)
{
char* cline = 0;
Elasticity* elp = dynamic_cast<Elasticity*>(myProblem);
if (!elp) return false;
if (!strncasecmp(keyWord,"GRAVITY",7))
{
double gx = atof(strtok(keyWord+7," "));
double gy = atof(strtok(NULL," "));
elp->setGravity(gx,gy,0.0);
if (myPid == 0)
std::cout <<"\nGravitation vector: "
<< gx <<" "<< gy << std::endl;
}
else if (!strncasecmp(keyWord,"ISOTROPHIC",10))
{
int nmat = atoi(keyWord+10);
if (myPid == 0)
std::cout <<"\nNumber of isotrophic materials: "<< nmat << std::endl;
for (int i = 0; i < nmat && (cline = utl::readLine(is)); i++)
{
int code = atoi(strtok(cline," "));
double E = atof(strtok(NULL," "));
double nu = atof(strtok(NULL," "));
double rho = atof(strtok(NULL," "));
if (myPid == 0)
std::cout <<"\tMaterial code "<< code <<": "
<< E <<" "<< nu <<" "<< rho << std::endl;
if (code == 0 || i == 0)
elp->setMaterial(E,nu,rho);
for (unsigned int j = 0; j < myProps.size() && code > 0; j++)
if (myProps[j].pindx == (size_t)code &&
myProps[j].pcode == Property::UNDEFINED)
{
myProps[j].pindx = mVec.size();
myProps[j].pcode = Property::MATERIAL;
}
if (code > 0)
mVec.push_back(IsoMat(E,nu,rho));
}
}
else if (!strncasecmp(keyWord,"CONSTANT_PRESSURE",17))
{
int npres = atoi(keyWord+17);
if (myPid == 0)
std::cout <<"\nNumber of pressures: "<< npres << std::endl;
for (int i = 0; i < npres && (cline = utl::readLine(is)); i++)
{
int code = atoi(strtok(cline," "));
int pdir = atoi(strtok(NULL," "));
double p = atof(strtok(NULL," "));
if (myPid == 0)
std::cout <<"\tPressure code "<< code <<" direction "<< pdir
<<": "<< p << std::endl;
this->setPropertyType(code,Property::NEUMANN);
myTracs[code] = new PressureField(p,pdir);
}
}
else if (!strncasecmp(keyWord,"ANASOL",6))
{
cline = strtok(keyWord+6," ");
if (!strncasecmp(cline,"HOLE",4))
{
double a = atof(strtok(NULL," "));
double F0 = atof(strtok(NULL," "));
double nu = atof(strtok(NULL," "));
asol = new Hole(a,F0,nu);
std::cout <<"\nAnalytical solution: Hole a="<< a <<" F0="<< F0
<<" nu="<< nu << std::endl;
}
else if (!strncasecmp(cline,"LSHAPE",6))
{
double a = atof(strtok(NULL," "));
double F0 = atof(strtok(NULL," "));
double nu = atof(strtok(NULL," "));
asol = new Lshape(a,F0,nu);
std::cout <<"\nAnalytical solution: Lshape a="<< a <<" F0="<< F0
<<" nu="<< nu << std::endl;
}
else if (!strncasecmp(cline,"CANTS",5))
{
double L = atof(strtok(NULL," "));
double H = atof(strtok(NULL," "));
double F0 = atof(strtok(NULL," "));
asol = new CanTS(L,H,F0);
std::cout <<"\nAnalytical solution: CanTS L="<< L <<" H="<< H
<<" F0="<< F0 << std::endl;
}
else if (!strncasecmp(cline,"CANTM",5))
{
double H = atof(strtok(NULL," "));
double M0 = atof(strtok(NULL," "));
asol = new CanTM(H,M0);
std::cout <<"\nAnalytical solution: CanTM H="<< H
<<" M0="<< M0 << std::endl;
}
else if (!strncasecmp(cline,"CURVED",6))
{
double a = atof(strtok(NULL," "));
double b = atof(strtok(NULL," "));
double u0 = atof(strtok(NULL," "));
double E = atof(strtok(NULL," "));
asol = new CurvedBeam(u0,a,b,E);
std::cout <<"\nAnalytical solution: Curved Beam a="<< a <<" b="<< b
<<" u0="<< u0 <<" E="<< E << std::endl;
}
else
std::cerr <<" ** SIMLinEl2D::parse: Unknown analytical solution "
<< cline << std::endl;
// Define the analytical boundary traction field
int code = (cline = strtok(NULL," ")) ? atoi(cline) : 0;
if (code > 0 && asol)
{
std::cout <<"Pressure code "<< code <<": Analytical traction"<< std::endl;
this->setPropertyType(code,Property::NEUMANN);
myTracs[code] = new TractionField(*asol);
}
}
// The remaining keywords are retained for backward compatibility with the
// prototype version. They enable direct specification of properties onto
// the topological entities (blocks and faces) of the model.
else if (!strncasecmp(keyWord,"PRESSURE",8))
{
Property press;
press.pcode = Property::NEUMANN;
press.ldim = 1;
int npres = atoi(keyWord+8);
std::cout <<"\nNumber of pressures: "<< npres << std::endl;
for (int i = 0; i < npres && (cline = utl::readLine(is)); i++)
{
press.pindx = 1+i;
press.patch = atoi(strtok(cline," "));
int pid = this->getLocalPatchIndex(press.patch);
if (pid < 0) return false;
if (pid < 1) continue;
press.lindx = atoi(strtok(NULL," "));
if (press.lindx < 1 || press.lindx > 4)
{
std::cerr <<" *** SIMLinEl2D::parse: Invalid edge index "
<< (int)press.lindx << std::endl;
return false;
}
if (asol)
{
std::cout <<"\tTraction on P"<< press.patch
<<" E"<< (int)press.lindx << std::endl;
myTracs[1+i] = new TractionField(*asol);
}
else
{
int pdir = atoi(strtok(NULL," "));
double p = atof(strtok(NULL," "));
std::cout <<"\tPressure on P"<< press.patch
<<" E"<< (int)press.lindx <<" direction "<< pdir <<": ";
if ((cline = strtok(NULL," ")))
myTracs[1+i] = new PressureField(utl::parseRealFunc(cline,p),pdir);
else
{
std::cout << p;
myTracs[1+i] = new PressureField(p,pdir);
}
std::cout << std::endl;
}
press.patch = pid;
myProps.push_back(press);
}
}
else if (!strncasecmp(keyWord,"MATERIAL",8))
{
int nmat = atoi(keyWord+8);
std::cout <<"\nNumber of materials: "<< nmat << std::endl;
for (int i = 0; i < nmat && (cline = utl::readLine(is)); i++)
{
double E = atof(strtok(cline," "));
double nu = atof(strtok(NULL," "));
double rho = atof(strtok(NULL," "));
while ((cline = strtok(NULL," ")))
if (!strncasecmp(cline,"ALL",3))
{
std::cout <<"\tMaterial for all patches: "
<< E <<" "<< nu <<" "<< rho << std::endl;
elp->setMaterial(E,nu,rho);
}
else
{
int patch = atoi(cline);
int pid = this->getLocalPatchIndex(patch);
if (pid < 0) return false;
if (pid < 1) continue;
std::cout <<"\tMaterial for P"<< patch
<<": "<< E <<" "<< nu <<" "<< rho << std::endl;
myProps.push_back(Property(Property::MATERIAL,mVec.size(),pid,3));
mVec.push_back(IsoMat(E,nu,rho));
if (i == 0)
elp->setMaterial(E,nu,rho);
}
}
}
else
return this->SIM2D::parse(keyWord,is);
return true;
}
bool SIMLinEl2D::initMaterial (size_t propInd)
{
Elasticity* elp = dynamic_cast<Elasticity*>(myProblem);
if (!elp) return false;
if (propInd >= mVec.size()) propInd = mVec.size()-1;
elp->setMaterial(mVec[propInd].E,mVec[propInd].nu,mVec[propInd].rho);
return true;
}
bool SIMLinEl2D::initNeumann (size_t propInd)
{
TracFuncMap::const_iterator tit = myTracs.find(propInd);
if (tit == myTracs.end()) return false;
Elasticity* elp = dynamic_cast<Elasticity*>(myProblem);
if (!elp) return false;
elp->setTraction(tit->second);
return true;
}

74
Apps/SIMLinEl2D.h Normal file
View File

@@ -0,0 +1,74 @@
// $Id: SIMLinEl2D.h,v 1.9 2010-12-18 16:23:52 kmo Exp $
//==============================================================================
//!
//! \file SIMLinEl2D.h
//!
//! \date Feb 04 2010
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Solution driver for 2D NURBS-based linear elastic FEM analysis.
//!
//==============================================================================
#ifndef _SIM_LIN_EL_2D_H
#define _SIM_LIN_EL_2D_H
#include "SIM2D.h"
#include "SIMenums.h"
/*!
\brief Driver class for 2D isogeometric FEM analysis of elasticity problems.
\details The class incapsulates data and methods for solving linear elasticity
problems using NURBS-based finite elements. It reimplements the parse method
of the parent class.
*/
class SIMLinEl2D : public SIM2D
{
//! \brief Struct for storage of physical material property parameters.
struct IsoMat
{
double E; //!< Young's modulus
double nu; //!< Poisson's ratio
double rho; //!< Mass density
//! \brief Default constructor.
IsoMat() { E = nu = rho = 0.0; }
//! \brief Constructor initializing a material instance.
IsoMat(double Emod, double Poiss, double D) : E(Emod), nu(Poiss), rho(D) {}
};
typedef std::vector<IsoMat> MaterialVec; //!< Vector of material data sets
public:
//! \brief Default constructor.
//! \param[in] form Problem formulation option
//! \param[in] planeStress Plane stress/plane strain option
SIMLinEl2D(int form = SIM::LINEAR, bool planeStress = true);
//! \brief The destructor frees the dynamically allocated data.
virtual ~SIMLinEl2D();
protected:
//! \brief Parses a data section from the input stream.
//! \param[in] keyWord Keyword of current data section to read
//! \param is The file stream to read from
virtual bool parse(char* keyWord, std::istream& is);
//! \brief Initializes material properties for integration of interior terms.
//! \param[in] propInd Physical property index
virtual bool initMaterial(size_t propInd);
//! \brief Initializes for integration of Neumann terms for a given property.
//! \param[in] propInd Physical property index
virtual bool initNeumann(size_t propInd);
//! \brief Returns the analytical solution function, if any.
virtual TensorFunc* getAnaSol() const { return asol; }
private:
MaterialVec mVec; //!< Material data
TensorFunc* asol; //!< Analytical stress field
};
#endif

428
Apps/SIMLinEl3D.C Normal file
View File

@@ -0,0 +1,428 @@
// $Id: SIMLinEl3D.C,v 1.28 2011-02-08 09:06:02 kmo Exp $
//==============================================================================
//!
//! \file SIMLinEl3D.C
//!
//! \date Dec 08 2009
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Solution driver for 3D NURBS-based linear elastic FEM analysis.
//!
//==============================================================================
#include "SIMLinEl3D.h"
#include "LinearElasticity.h"
#include "NonlinearElasticity.h"
#include "AnalyticSolutions.h"
#include "Functions.h"
#include "Utilities.h"
#include "Vec3Oper.h"
#include "Property.h"
#include "Tensor.h"
#include <string.h>
/*!
\brief Local coordinate system for a cylinder along global z-axis.
*/
class CylinderCS : public LocalSystem
{
public:
//! \brief Constructor printing a message making user aware of its presense.
CylinderCS()
{
std::cout <<"\nLocal coordinate system: Cylindric"<< std::endl;
}
//! \brief Computes the global-to-local transformation at the point \a X.
virtual const Tensor& getTmat(const Vec3& X) const
{
static Tensor T(3);
double r = hypot(X.x,X.y);
T(1,1) = X.x/r;
T(1,2) = X.y/r;
T(2,1) = -T(1,2);
T(2,2) = T(1,1);
T(3,3) = 1.0;
return T;
}
};
#ifdef PRINT_CS
#include <fstream>
#endif
/*!
\brief Local coordinate system for a cylinder along global z-axis,
closed by a spherical cap.
*/
class CylinderSphereCS : public LocalSystem
{
public:
//! \brief Constructor printing a message making user aware of its presense.
CylinderSphereCS(double H = 0.0) : h(H)
{
std::cout <<"\nLocal coordinate system: Cylindric with Spherical cap, h="
<< h << std::endl;
#ifdef PRINT_CS
sn.open("nodes.dat");
se.open("elements.dat");
s1.open("v1.dat");
s2.open("v2.dat");
s3.open("v3.dat");
sn <<"\n*NODES 4\n";
se <<"\n*ELEMENTS 4\n%NODES #4\n"
<<"%NO_ID\n%MAP_NODE_INDICES\n%PART_ID 4\n%POINTS\n";
s1 <<"\n*RESULTS 31\n%NO_ID\n%DIMENSION 3\n%PER_NODE #4\n";
s2 <<"\n*RESULTS 32\n%NO_ID\n%DIMENSION 3\n%PER_NODE #4\n";
s3 <<"\n*RESULTS 33\n%NO_ID\n%DIMENSION 3\n%PER_NODE #4\n";
}
virtual ~CylinderSphereCS()
{
std::cout <<"Finalizing VTF-output of local coordinate systems"<< std::endl;
s1 <<"\n*GLVIEWVECTOR 2\n%NAME \"v1\"\n%STEP 1\n31\n";
s2 <<"\n*GLVIEWVECTOR 3\n%NAME \"v2\"\n%STEP 1\n32\n";
s3 <<"\n*GLVIEWVECTOR 4\n%NAME \"v3\"\n%STEP 1\n33\n";
#endif
}
//! \brief Computes the global-to-local transformation at the point \a X.
virtual const Tensor& getTmat(const Vec3& X) const
{
static Tensor T(3);
#ifdef PRINT_CS
sn << X <<'\n';
static int iel = 0;
se << ++iel <<'\n';
#endif
if (patch == 1) // Cylindric system {-z,theta,r}
{
T.zero();
double r = hypot(X.x,X.y);
T(1,3) = -1.0;
T(2,1) = -X.y/r;
T(2,2) = X.x/r;
T(3,1) = T(2,2);
T(3,2) = -T(2,1);
#ifdef PRINT_CS
s1 <<"0 0 -1\n";
s2 << T(2,1) <<" "<< T(2,2) <<" 0\n";
s3 << T(3,1) <<" "<< T(3,2) <<" 0\n";
#endif
}
else // Spherical system {phi,theta,r}
{
Vec3 v3(X.x,X.y,X.z-h);
v3 /= v3.length();
double theta = atan2(X.y,X.x);
double phi = acos(v3.z);
Vec3 v1(cos(theta)*cos(phi),sin(theta)*cos(phi),-sin(phi));
Vec3 v2(v3,v1);
for (int i = 1; i <= 3; i++)
{
T(1,i) = v1[i-1];
T(2,i) = v2[i-1];
T(3,i) = v3[i-1];
}
#ifdef PRINT_CS
s1 << v1 <<'\n';
s2 << v2 <<'\n';
s3 << v3 <<'\n';
#endif
}
return T;
}
private:
double h; //!< Height above global origin of the centre of the sphere
#ifdef PRINT_CS
mutable std::ofstream sn; //!< VTF output stream for CS nodes
mutable std::ofstream se; //!< VTF output stream for CS point elements
mutable std::ofstream s1; //!< VTF output stream for vector v1 of local CS
mutable std::ofstream s2; //!< VTF output stream for vector v2 of local CS
mutable std::ofstream s3; //!< VTF output stream for vector v3 of local CS
#endif
};
SIMLinEl3D::SIMLinEl3D (bool checkRHS, int form) : SIM3D(checkRHS)
{
switch (form)
{
case SIM::LINEAR:
myProblem = new LinearElasticity();
break;
case SIM::NONLINEAR:
myProblem = new NonlinearElasticity();
break;
}
asol = 0;
}
SIMLinEl3D::~SIMLinEl3D ()
{
if (asol) delete asol;
}
bool SIMLinEl3D::parse (char* keyWord, std::istream& is)
{
char* cline = 0;
int nConstPress = 0;
int nLinearPress = 0;
Elasticity* elp = dynamic_cast<Elasticity*>(myProblem);
if (!elp) return false;
if (!strncasecmp(keyWord,"GRAVITY",7))
{
double gx = atof(strtok(keyWord+7," "));
double gy = atof(strtok(NULL," "));
double gz = atof(strtok(NULL," "));
elp->setGravity(gx,gy,gz);
if (myPid == 0)
std::cout <<"\nGravitation vector: "
<< gx <<" "<< gy <<" "<< gz << std::endl;
}
else if (!strncasecmp(keyWord,"ISOTROPHIC",10))
{
int nmat = atoi(keyWord+10);
if (myPid == 0)
std::cout <<"\nNumber of isotrophic materials: "<< nmat << std::endl;
for (int i = 0; i < nmat && (cline = utl::readLine(is)); i++)
{
int code = atoi(strtok(cline," "));
double E = atof(strtok(NULL," "));
double nu = atof(strtok(NULL," "));
double rho = atof(strtok(NULL," "));
if (myPid == 0)
std::cout <<"\tMaterial code "<< code <<": "
<< E <<" "<< nu <<" "<< rho << std::endl;
if (code == 0 || i == 0)
elp->setMaterial(E,nu,rho);
for (unsigned int j = 0; j < myProps.size() && code > 0; j++)
if (myProps[j].pindx == code && myProps[j].pcode == Property::UNDEFINED)
{
myProps[j].pindx = mVec.size();
myProps[j].pcode = Property::MATERIAL;
}
if (code > 0)
mVec.push_back(IsoMat(E,nu,rho));
}
}
else if (!strncasecmp(keyWord,"CONSTANT_PRESSURE",17))
nConstPress = atoi(keyWord+17);
else if (!strncasecmp(keyWord,"LINEAR_PRESSURE",15))
nLinearPress = atoi(keyWord+15);
else if (!strncasecmp(keyWord,"ANASOL",6))
{
cline = strtok(keyWord+6," ");
if (!strncasecmp(cline,"HOLE",4))
{
double a = atof(strtok(NULL," "));
double F0 = atof(strtok(NULL," "));
double nu = atof(strtok(NULL," "));
asol = new Hole(a,F0,nu,true);
std::cout <<"\nAnalytical solution: Hole a="<< a <<" F0="<< F0
<<" nu="<< nu << std::endl;
}
else if (!strncasecmp(cline,"LSHAPE",6))
{
double a = atof(strtok(NULL," "));
double F0 = atof(strtok(NULL," "));
double nu = atof(strtok(NULL," "));
asol = new Lshape(a,F0,nu,true);
std::cout <<"\nAnalytical solution: Lshape a="<< a <<" F0="<< F0
<<" nu="<< nu << std::endl;
}
else if (!strncasecmp(cline,"CANTS",5))
{
double L = atof(strtok(NULL," "));
double H = atof(strtok(NULL," "));
double F0 = atof(strtok(NULL," "));
asol = new CanTS(L,H,F0,true);
std::cout <<"\nAnalytical solution: CanTS L="<< L <<" H="<< H
<<" F0="<< F0 << std::endl;
}
else
std::cerr <<" ** SIMLinEl3D::parse: Unknown analytical solution "
<< cline << std::endl;
// Define the analytical boundary traction field
int code = (cline = strtok(NULL," ")) ? atoi(cline) : 0;
if (code > 0 && asol)
{
std::cout <<"Pressure code "<< code <<": Analytical traction"<< std::endl;
this->setPropertyType(code,Property::NEUMANN);
myTracs[code] = new TractionField(*asol);
}
}
// The remaining keywords are retained for backward compatibility with the
// prototype version. They enable direct specification of properties onto
// the topological entities (blocks and faces) of the model.
else if (!strncasecmp(keyWord,"PRESSURE",8))
{
Property press;
press.pcode = Property::NEUMANN;
press.ldim = 2;
int npres = atoi(keyWord+8);
std::cout <<"\nNumber of pressures: "<< npres << std::endl;
for (int i = 0; i < npres && (cline = utl::readLine(is)); i++)
{
press.pindx = 1+i;
press.patch = atoi(strtok(cline," "));
int pid = this->getLocalPatchIndex(press.patch);
if (pid < 0) return false;
if (pid < 1) continue;
press.lindx = atoi(strtok(NULL," "));
if (press.lindx < 1 || press.lindx > 6)
{
std::cerr <<" *** SIMLinEl3D::parse: Invalid face index "
<< (int)press.lindx << std::endl;
return false;
}
if (asol)
{
std::cout <<"\tTraction on P"<< press.patch
<<" F"<< (int)press.lindx << std::endl;
myTracs[1+i] = new TractionField(*asol);
}
else
{
int pdir = atoi(strtok(NULL," "));
double p = atof(strtok(NULL," "));
std::cout <<"\tPressure on P"<< press.patch
<<" F"<< (int)press.lindx <<" direction "<< pdir <<": ";
if (cline = strtok(NULL," "))
myTracs[1+i] = new PressureField(utl::parseRealFunc(cline,p),pdir);
else
{
std::cout << p;
myTracs[1+i] = new PressureField(p,pdir);
}
std::cout << std::endl;
}
press.patch = pid;
myProps.push_back(press);
}
}
else if (!strncasecmp(keyWord,"MATERIAL",8))
{
int nmat = atoi(keyWord+8);
std::cout <<"\nNumber of materials: "<< nmat << std::endl;
for (int i = 0; i < nmat && (cline = utl::readLine(is)); i++)
{
double E = atof(strtok(cline," "));
double nu = atof(strtok(NULL," "));
double rho = atof(strtok(NULL," "));
while (cline = strtok(NULL," "))
if (!strncasecmp(cline,"ALL",3))
{
std::cout <<"\tMaterial for all patches: "
<< E <<" "<< nu <<" "<< rho << std::endl;
elp->setMaterial(E,nu,rho);
}
else
{
int patch = atoi(cline);
int pid = this->getLocalPatchIndex(patch);
if (pid < 0) return false;
if (pid < 1) continue;
std::cout <<"\tMaterial for P"<< patch
<<": "<< E <<" "<< nu <<" "<< rho << std::endl;
myProps.push_back(Property(Property::MATERIAL,mVec.size(),pid,3));
mVec.push_back(IsoMat(E,nu,rho));
if (i == 0)
elp->setMaterial(E,nu,rho);
}
}
}
else if (!strncasecmp(keyWord,"LOCAL_SYSTEM",12))
{
size_t i = 12;
while (i < strlen(keyWord) && isspace(keyWord[i])) i++;
if (!strncasecmp(keyWord+i,"CYLINDRICZ",10))
elp->setLocalSystem(new CylinderCS);
else if (!strncasecmp(keyWord+i,"CYLINDER+SPHERE",15))
elp->setLocalSystem(new CylinderSphereCS(atof(keyWord+i+15)));
else
std::cerr <<" *** SIMLinEl3D::parse: Unsupported coordinate system: "
<< keyWord+i << std::endl;
}
else
return this->SIM3D::parse(keyWord,is);
int npres = nConstPress + nLinearPress;
if (npres > 0)
{
if (myPid == 0)
std::cout <<"\nNumber of pressures: "<< npres << std::endl;
for (int i = 0; i < npres && (cline = utl::readLine(is)); i++)
{
int code = atoi(strtok(cline," "));
int pdir = atoi(strtok(NULL," "));
double p = atof(strtok(NULL," "));
if (myPid == 0)
std::cout <<"\tPressure code "<< code <<" direction "<< pdir
<<": "<< p << std::endl;
this->setPropertyType(code,Property::NEUMANN);
if (nLinearPress)
{
RealFunc* pfl = new ConstTimeFunc(new LinearFunc(p));
myTracs[code] = new PressureField(pfl,pdir);
}
else
myTracs[code] = new PressureField(p,pdir);
}
}
return true;
}
bool SIMLinEl3D::initMaterial (size_t propInd)
{
Elasticity* elp = dynamic_cast<Elasticity*>(myProblem);
if (!elp) return false;
if (propInd >= mVec.size()) propInd = mVec.size()-1;
elp->setMaterial(mVec[propInd].E,mVec[propInd].nu,mVec[propInd].rho);
return true;
}
bool SIMLinEl3D::initNeumann (size_t propInd)
{
TracFuncMap::const_iterator tit = myTracs.find(propInd);
if (tit == myTracs.end()) return false;
Elasticity* elp = dynamic_cast<Elasticity*>(myProblem);
if (!elp) return false;
elp->setTraction(tit->second);
return true;
}

74
Apps/SIMLinEl3D.h Normal file
View File

@@ -0,0 +1,74 @@
// $Id: SIMLinEl3D.h,v 1.11 2010-12-18 16:23:52 kmo Exp $
//==============================================================================
//!
//! \file SIMLinEl3D.h
//!
//! \date Dec 08 2009
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Solution driver for 3D NURBS-based linear elastic FEM analysis.
//!
//==============================================================================
#ifndef _SIM_LIN_EL_3D_H
#define _SIM_LIN_EL_3D_H
#include "SIM3D.h"
#include "SIMenums.h"
/*!
\brief Driver class for 3D isogeometric FEM analysis of elasticity problems.
\details The class incapsulates data and methods for solving linear elasticity
problems using NURBS-based finite elements. It reimplements the parse method
of the parent class.
*/
class SIMLinEl3D : public SIM3D
{
//! \brief Struct for storage of physical material property parameters.
struct IsoMat
{
double E; //!< Young's modulus
double nu; //!< Poisson's ratio
double rho; //!< Mass density
//! \brief Default constructor.
IsoMat() { E = nu = rho = 0.0; }
//! \brief Constructor initializing a material instance.
IsoMat(double Emod, double Poiss, double D) : E(Emod), nu(Poiss), rho(D) {}
};
typedef std::vector<IsoMat> MaterialVec; //!< Vector of material data sets
public:
//! \brief Default constructor.
//! \param[in] checkRHS If \e true, ensure the model is in a right-hand system
//! \param[in] form Problem formulation option
SIMLinEl3D(bool checkRHS = false, int form = SIM::LINEAR);
//! \brief The destructor frees the dynamically allocated data.
virtual ~SIMLinEl3D();
protected:
//! \brief Parses a data section from the input stream.
//! \param[in] keyWord Keyword of current data section to read
//! \param is The file stream to read from
virtual bool parse(char* keyWord, std::istream& is);
//! \brief Initializes material properties for integration of interior terms.
//! \param[in] propInd Physical property index
virtual bool initMaterial(size_t propInd);
//! \brief Initializes for integration of Neumann terms for a given property.
//! \param[in] propInd Physical property index
virtual bool initNeumann(size_t propInd);
//! \brief Returns the analytical solution function, if any.
virtual TensorFunc* getAnaSol() const { return asol; }
private:
MaterialVec mVec; //!< Material data
TensorFunc* asol; //!< Analytical stress field
};
#endif

386
Apps/main_LinEl3D.C Normal file
View File

@@ -0,0 +1,386 @@
// $Id: main_LinEl3D.C,v 1.27 2011-02-08 09:32:18 kmo Exp $
//==============================================================================
//!
//! \file main_LinEl3D.C
//!
//! \date Jan 28 2009
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Main program for the isogeometric linear elasticity solver.
//!
//==============================================================================
#include "SIMLinEl3D.h"
#include "SIMLinEl2D.h"
#include "LinAlgInit.h"
#include "Profiler.h"
#include <fstream>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
/*!
\brief Main program for the NURBS-based isogeometric linear elasticity solver.
The input to the program is specified through the following
command-line arguments. The arguments may be given in arbitrary order.
\arg \a input-file : Input file with model definition
\arg -dense : Use the dense LAPACK matrix equation solver
\arg -spr : Use the SPR direct equation solver
\arg -superlu : Use the sparse SuperLU equation solver
\arg -samg : Use the sparse algebraic multi-grid equation solver
\arg -petsc : Use equation solver from PETSc library
\arg -nGauss \a n : Number of Gauss points over a knot-span in each direction
\arg -vtf \a format : VTF-file format (-1=NONE, 0=ASCII, 1=BINARY)
\arg -nviz \a nviz : Number of visualization points over each knot-span
\arg -nu \a nu : Number of visualization points per knot-span in u-direction
\arg -nv \a nv : Number of visualization points per knot-span in v-direction
\arg -nw \a nw : Number of visualization points per knot-span in w-direction
\arg -dumpASC : Dump model and solution to ASCII files for external processing
\arg -ignore \a p1, \a p2, ... : Ignore these patches in the analysis
\arg -eig \a iop : Eigenproblem solver to use (1...6)
\arg -nev \a nev : Number of eigenvalues to compute
\arg -ncv \a ncv : Number of Arnoldi vectors to use in the eigenvalue analysis
\arg -shift \a shf : Shift value to use in the eigenproblem solver
\arg -free : Ignore all boundary conditions (use in free vibration analysis)
\arg -check : Data check only, read model and output to VTF (no solution)
\arg -checkRHS : Check that the patches are modelled in a right-hand system
\arg -vizRHS : Save the right-hand-side load vector on the VTF-file
\arg -fixDup : Resolve co-located nodes by merging them into a single node
\arg -2D : Use two-parametric simulation driver (plane stress)
\arg -2Dpstrain : Use two-parametric simulation driver (plane strain)
\arg -TL : Use the nonlinear total-lagrangian elasticity integrand
\arg -lag : Use Lagrangian basis functions instead of splines/NURBS
\arg -spec : Use Spectral basis functions instead of splines/NURBS
*/
int main (int argc, char** argv)
{
Profiler prof(argv[0]);
utl::profiler->start("Initialization");
SIM::Formulation form = SIM::LINEAR;
SystemMatrix::Type solver = SystemMatrix::SPARSE;
int nGauss = 4;
int format = -1;
int n[3] = { 2, 2, 2 };
std::vector<int> ignoredPatches;
int iop = 0;
int nev = 10;
int ncv = 20;
double shf = 0.0;
bool checkRHS = false;
bool vizRHS = false;
bool fixDup = false;
bool dumpASCII = false;
char twoD = false;
char* infile = 0;
LinAlgInit linalg(argc,argv);
for (int i = 1; i < argc; i++)
if (!strcmp(argv[i],"-dense"))
solver = SystemMatrix::DENSE;
else if (!strcmp(argv[i],"-spr"))
solver = SystemMatrix::SPR;
else if (!strncmp(argv[i],"-superlu",8))
{
solver = SystemMatrix::SPARSE;
if (isdigit(argv[i][8]))
SIMbase::num_threads_SLU = atoi(argv[i]+8);
}
else if (!strcmp(argv[i],"-samg"))
solver = SystemMatrix::SAMG;
else if (!strcmp(argv[i],"-petsc"))
solver = SystemMatrix::PETSC;
else if (!strcmp(argv[i],"-nGauss") && i < argc-1)
nGauss = atoi(argv[++i]);
else if (!strcmp(argv[i],"-vtf") && i < argc-1)
format = atoi(argv[++i]);
else if (!strcmp(argv[i],"-nviz") && i < argc-1)
n[0] = n[1] = n[2] = atoi(argv[++i]);
else if (!strcmp(argv[i],"-nu") && i < argc-1)
n[0] = atoi(argv[++i]);
else if (!strcmp(argv[i],"-nv") && i < argc-1)
n[1] = atoi(argv[++i]);
else if (!strcmp(argv[i],"-nw") && i < argc-1)
n[2] = atoi(argv[++i]);
else if (!strcmp(argv[i],"-dumpASC"))
dumpASCII = true;
else if (!strcmp(argv[i],"-ignore"))
while (i < argc-1 && isdigit(argv[i+1][0]))
{
char* endp = 0;
int endVal = 0;
ignoredPatches.push_back(strtol(argv[++i],&endp,10));
if (endp && *endp == ':')
endVal = strtol(endp+1,&endp,10);
while (ignoredPatches.back() < endVal)
ignoredPatches.push_back(ignoredPatches.back()+1);
}
else if (!strcmp(argv[i],"-eig") && i < argc-1)
iop = atoi(argv[++i]);
else if (!strcmp(argv[i],"-nev") && i < argc-1)
nev = atoi(argv[++i]);
else if (!strcmp(argv[i],"-ncv") && i < argc-1)
ncv = atoi(argv[++i]);
else if (!strcmp(argv[i],"-shift") && i < argc-1)
shf = atof(argv[++i]);
else if (!strcmp(argv[i],"-free"))
SIMbase::ignoreDirichlet = true;
else if (!strcmp(argv[i],"-check"))
iop = 100;
else if (!strcmp(argv[i],"-checkRHS"))
checkRHS = true;
else if (!strcmp(argv[i],"-vizRHS"))
vizRHS = true;
else if (!strcmp(argv[i],"-fixDup"))
fixDup = true;
else if (!strncmp(argv[i],"-2D",3))
twoD = strcmp(argv[i],"-2Dpstrain") ? 1 : 2;
else if (!strcmp(argv[i],"-TL"))
form = SIM::NONLINEAR;
else if (!strncmp(argv[i],"-lag",4))
SIMbase::discretization = SIMbase::Lagrange;
else if (!strncmp(argv[i],"-spec",5))
SIMbase::discretization = SIMbase::Spectral;
else if (!infile)
infile = argv[i];
else
std::cerr <<" ** Unknown option ignored: "<< argv[i] << std::endl;
if (!infile)
{
std::cout <<"usage: "<< argv[0]
<<" <inputfile> [-dense|-spr|-superlu<nt>|-samg|-petsc]\n "
<<" [-free] [-lag] [-spec] [-2D[pstrain]] [-TL] [-nGauss <n>]\n"
<<" [-vtf <format>] [-nviz <nviz>]"
<<" [-nu <nu>] [-nv <nv>] [-nw <nw>]\n"
<<" [-eig <iop>] [-nev <nev>] [-ncv <ncv] [-shift <shf>]\n"
<<" [-ignore <p1> <p2> ...] [-fixDup]"
<<" [-checkRHS] [-check] [-dumpASC]\n";
return 0;
}
if (twoD) n[2] = 1;
// Load vector visualization is not available when using additional viz-points
if (n[0] > 2 || n[1] > 2 || n[2] > 2) vizRHS = false;
// Boundary conditions can be ignored only in generalized eigenvalue analysis
if (iop != 4 && iop != 6) SIMbase::ignoreDirichlet = false;
if (linalg.myPid == 0)
{
std::cout <<"\n >>> Spline FEM Linear Elasticity solver <<<"
<<"\n ===========================================\n"
<<"\nInput file: "<< infile
<<"\nEquation solver: "<< solver
<<"\nNumber of Gauss points: "<< nGauss;
if (format >= 0)
{
std::cout <<"\nVTF file format: "<< (format ? "BINARY":"ASCII")
<<"\nNumber of visualization points: "
<< n[0] <<" "<< n[1];
if (!twoD) std::cout <<" "<< n[2];
}
if (iop > 0 && iop < 100)
std::cout <<"\nEigenproblem solver: "<< iop
<<"\nNumber of eigenvalues: "<< nev
<<"\nNumber of Arnoldi vectors: "<< ncv
<<"\nShift value: "<< shf;
if (form == SIM::NONLINEAR)
std::cout <<"\nGreen-Lagrange strains are used";
if (SIMbase::discretization == SIMbase::Lagrange)
std::cout <<"\nLagrangian basis functions are used";
else if (SIMbase::discretization == SIMbase::Spectral)
std::cout <<"\nSpectral basis functions are used";
if (SIMbase::ignoreDirichlet)
std::cout <<"\nSpecified boundary conditions are ignored";
if (fixDup)
std::cout <<"\nCo-located nodes will be merged";
if (checkRHS)
std::cout <<"\nCheck that each patch has a right-hand coordinate system";
if (!ignoredPatches.empty())
{
std::cout <<"\nIgnored patches:";
for (size_t i = 0; i < ignoredPatches.size(); i++)
std::cout <<" "<< ignoredPatches[i];
}
std::cout << std::endl;
}
utl::profiler->stop("Initialization");
utl::profiler->start("Model input");
// Read in model definitions and establish the FE data structures
SIMbase* model;
if (twoD)
model = new SIMLinEl2D(form,twoD==1);
else
model = new SIMLinEl3D(checkRHS,form);
if (!model->read(infile))
return 1;
utl::profiler->stop("Model input");
if (linalg.myPid == 0)
model->printProblem(std::cout);
if (!model->preprocess(ignoredPatches,fixDup))
return 1;
model->setQuadratureRule(nGauss);
Matrix eNorm;
Vector gNorm, load;
Vectors displ(1);
std::vector<Mode> modes;
std::vector<Mode>::const_iterator it;
switch (iop) {
case 0:
case 5:
// Static solution: Assemble [Km] and {R}
model->setMode(SIM::STATIC);
model->initSystem(solver,1,1);
model->setAssociatedRHS(0,0);
if (!model->assembleSystem())
return 2;
else if (vizRHS)
model->extractLoadVec(load);
// Solve the linear system of equations
if (!model->solveSystem(displ.front(),1))
return 3;
// Evaluate solution norms
model->setMode(SIM::RECOVERY);
if (!model->solutionNorms(displ,eNorm,gNorm))
return 4;
if (linalg.myPid == 0)
{
std::cout <<"Energy norm |u^h| = a(u^h,u^h)^0.5 : "<< gNorm(1);
std::cout <<"\nExternal energy ((f,u^h)+(t,u^h)^0.5 : "<< gNorm(2);
if (gNorm.size() > 2)
std::cout <<"\nExact norm |u| = a(u,u)^0.5 : "<< gNorm(3);
if (gNorm.size() > 3)
std::cout <<"\nExact error a(e,e)^0.5, e=u-u^h : "<< gNorm(4)
<<"\nExact relative error (%) : "<< gNorm(4)/gNorm(3)*100.0;
std::cout << std::endl;
}
if (iop == 0) break;
// Linearized buckling: Assemble [Km] and [Kg]
model->setMode(SIM::BUCKLING);
model->initSystem(solver,2,0);
if (!model->assembleSystem(displ))
return 5;
// Solve the generalized eigenvalue problem
if (!model->systemModes(modes,nev,ncv,iop,shf))
return 6;
break;
case 1:
case 2:
// Assemble and solve the regular eigenvalue problem
model->setMode(SIM::STIFF_ONLY);
model->initSystem(solver,1,0);
if (!model->assembleSystem())
return 5;
if (!model->systemModes(modes,nev,ncv,iop,shf))
return 6;
break;
case 100:
break; // Model check
default:
// Free vibration: Assemble [Km] and [M]
model->setMode(SIM::VIBRATION);
model->initSystem(solver,2,0);
if (!model->assembleSystem())
return 5;
if (!model->systemModes(modes,nev,ncv,iop,shf))
return 6;
}
utl::profiler->start("Postprocessing");
if (format >= 0)
{
// Write VTF-file with model geometry
int iStep = 1, nBlock = 0;
if (!model->writeGlv(infile,n,format))
return 7;
// Write boundary tractions, if any
if (!model->writeGlvT(iStep,nBlock))
return 8;
// Write Dirichlet boundary conditions
if (!model->writeGlvBC(n,nBlock))
return 8;
// Write load vector to VTF-file
if (!model->writeGlvV(load,"Load vector",n,iStep,nBlock))
return 9;
// Write solution fields to VTF-file
if (!model->writeGlvS(displ.front(),n,iStep,nBlock))
return 10;
// Write eigenmodes
for (it = modes.begin(); it != modes.end(); it++)
if (!model->writeGlvM(*it, iop==3 || iop==4 || iop==6, n, nBlock))
return 11;
// Write element norms (when no additional visualization points are used)
if (n[0] == 2 && n[1] == 2 && n[2] <= 2)
if (!model->writeGlvN(eNorm,iStep,nBlock))
return 12;
model->closeGlv();
}
if (dumpASCII)
{
// Write (refined) model to g2-file
std::ofstream osg(strcat(strtok(infile,"."),".g2"));
std::cout <<"\nWriting updated g2-file "<< infile << std::endl;
model->dumpGeometry(osg);
if (!displ.front().empty())
{
// Write solution (control point values) to ASCII files
std::ofstream osd(strcat(strtok(infile,"."),".dis"));
std::cout <<"\nWriting deformation to file "<< infile << std::endl;
model->dumpPrimSol(displ.front(),osd,false);
std::ofstream oss(strcat(strtok(infile,"."),".sol"));
std::cout <<"\nWriting solution to file "<< infile << std::endl;
model->dumpSolution(displ.front(),oss);
}
if (!modes.empty())
{
// Write eigenvectors to ASCII files
std::ofstream ose(strcat(strtok(infile,"."),".eig"));
std::cout <<"\nWriting eigenvectors to file "<< infile << std::endl;
for (it = modes.begin(); it != modes.end(); it++)
{
ose <<"# Eigenvector_"<< it->eigNo <<" Eigenvalue="<< it->eigVal <<"\n";
model->dumpPrimSol(it->eigVec,ose,false);
}
}
}
utl::profiler->stop("Postprocessing");
delete model;
return 0;
}

251
Apps/main_LinElStep.C Normal file
View File

@@ -0,0 +1,251 @@
// $Id: main_LinElStep.C,v 1.7 2011-02-05 18:32:25 kmo Exp $
//==============================================================================
//!
//! \file main_LinElStep.C
//!
//! \date Jun 1 2010
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Main program for the isogeometric linear elasticity solver.
//!
//==============================================================================
#include "SIMLinEl3D.h"
#include "LinAlgInit.h"
#include "TimeDomain.h"
#include "Profiler.h"
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
/*!
\brief Main program for the NURBS-based isogeometric linear elasticity solver.
This program is basicly the same as main_LinEl3D.C, except that this one runs
in a loop over a specified number of load steps to compute the linear solution
for each step, when the load function is assumed to vary with time.
This version is for 3D patches only, and no eigenvalue solution is included.
The input to the program is specified through the following
command-line arguments. The arguments may be given in arbitrary order.
\arg \a input-file : Input file with model definition
\arg -dense : Use the dense LAPACK matrix equation solver
\arg -spr : Use the SPR direct equation solver
\arg -superlu : Use the sparse superLU equation solver
\arg -samg : Use the sparse algebraic multi-grid equation solver
\arg -petsc : Use equation solver from PETSc library
\arg -nGauss \a n : Number of Gauss points over a knot-span in each direction
\arg -vtf \a format : VTF-file format (0=ASCII, 1=BINARY)
\arg -nviz \a nviz : Number of visualization points over each knot-span
\arg -nu \a nu : Number of visualization points per knot-span in u-direction
\arg -nv \a nv : Number of visualization points per knot-span in v-direction
\arg -nw \a nw : Number of visualization points per knot-span in w-direction
\arg -ignore \a p1, \a p2, ... : Ignore these patches in the analysis
\arg -checkRHS : Check that the patches are modelled in a right-hand system
\arg -vizRHS : Save the right-hand-side load vector on the VTF-file
\arg -fixDup : Resolve co-located nodes by merging them into a single node
\arg -lag : Use Lagrangian basis functions instead of splines/NURBS
\arg -spec : Use Spectral basis functions instead of splines/NURBS
\arg -nstep \a n : Number of time steps
\arg -dt \a dt : Time increment size
*/
int main (int argc, char** argv)
{
Profiler prof(argv[0]);
utl::profiler->start("Initialization");
SystemMatrix::Type solver = SystemMatrix::SPARSE;
int nGauss = 4;
int format = 0;
int n[3] = { 2, 2, 2 };
std::vector<int> ignoredPatches;
bool checkRHS = false;
bool vizRHS = false;
bool fixDup = false;
char* infile = 0;
TimeDomain time;
int nStep = 1;
LinAlgInit linalg(argc,argv);
for (int i = 1; i < argc; i++)
if (!strcmp(argv[i],"-dense"))
solver = SystemMatrix::DENSE;
else if (!strcmp(argv[i],"-spr"))
solver = SystemMatrix::SPR;
else if (!strcmp(argv[i],"-superlu"))
solver = SystemMatrix::SPARSE;
else if (!strcmp(argv[i],"-samg"))
solver = SystemMatrix::SAMG;
else if (!strcmp(argv[i],"-petsc"))
solver = SystemMatrix::PETSC;
else if (!strcmp(argv[i],"-nGauss") && i < argc-1)
nGauss = atoi(argv[++i]);
else if (!strcmp(argv[i],"-vtf") && i < argc-1)
format = atoi(argv[++i]);
else if (!strcmp(argv[i],"-nviz") && i < argc-1)
n[0] = n[1] = n[2] = atoi(argv[++i]);
else if (!strcmp(argv[i],"-nu") && i < argc-1)
n[0] = atoi(argv[++i]);
else if (!strcmp(argv[i],"-nv") && i < argc-1)
n[1] = atoi(argv[++i]);
else if (!strcmp(argv[i],"-nw") && i < argc-1)
n[2] = atoi(argv[++i]);
else if (!strcmp(argv[i],"-ignore"))
while (i < argc-1 && isdigit(argv[i+1][0]))
{
char* endp = 0;
int endVal = 0;
ignoredPatches.push_back(strtol(argv[++i],&endp,10));
if (endp && *endp == ':')
endVal = strtol(endp+1,&endp,10);
while (ignoredPatches.back() < endVal)
ignoredPatches.push_back(ignoredPatches.back()+1);
}
else if (!strcmp(argv[i],"-checkRHS"))
checkRHS = true;
else if (!strcmp(argv[i],"-vizRHS"))
vizRHS = true;
else if (!strcmp(argv[i],"-fixDup"))
fixDup = true;
else if (!strncmp(argv[i],"-lag",4))
SIMbase::discretization = SIMbase::Lagrange;
else if (!strncmp(argv[i],"-spec",5))
SIMbase::discretization = SIMbase::Spectral;
else if (!strcmp(argv[i],"-nstep") && i < argc-1)
nStep = atoi(argv[++i]);
else if (!strcmp(argv[i],"-dt") && i < argc-1)
time.dt = atof(argv[++i]);
else if (!infile)
infile = argv[i];
else
std::cerr <<" ** Unknown option ignored: "<< argv[i] << std::endl;
if (!infile)
{
std::cout <<"usage: "<< argv[0]
<<" <inputfile> [-dense|-spr|-superlu|-samg|-petsc]\n"
<<" [-lag] [-spec] [-nGauss <n>] [-nstep <n>] [-dt <dt>]\n"
<<" [-vtf <format>] [-nviz <nviz>]"
<<" [-nu <nu>] [-nv <nv>] [-nw <nw>]\n"
<<" [-ignore <p1> <p2> ...] [-fixDup] [-checkRHS]\n";
return 0;
}
// Load vector visualization is not available when using additional viz-points
if (n[0] > 2 || n[1] > 2 || n[2] > 2) vizRHS = false;
std::cout <<"\n >>> Spline FEM Linear Elasticity solver <<<"
<<"\n ===========================================\n"
<<"\nInput file: "<< infile
<<"\nEquation solver: "<< solver
<<"\nNumber of Gauss points: "<< nGauss
<<"\nVTF file format: "<< (format ? "BINARY":"ASCII")
<<"\nNumber of visualization points: "
<< n[0] <<" "<< n[1] <<" "<< n[2];
if (SIMbase::discretization == SIMbase::Lagrange)
std::cout <<"\nLagrangian basis functions are used";
else if (SIMbase::discretization == SIMbase::Spectral)
std::cout <<"\nSpectral basis functions are used";
if (fixDup)
std::cout <<"\nCo-located nodes will be merged";
if (checkRHS)
std::cout <<"\nCheck that each patch has a right-hand coordinate system";
if (!ignoredPatches.empty())
{
std::cout <<"\nIgnored patches:";
for (size_t i = 0; i < ignoredPatches.size(); i++)
std::cout <<" "<< ignoredPatches[i];
}
std::cout << std::endl;
utl::profiler->stop("Initialization");
utl::profiler->start("Model input");
// Read in model definitions and establish the FE structures
SIMbase* model = new SIMLinEl3D(checkRHS);
if (!model->read(infile) || !model->preprocess(ignoredPatches,fixDup))
return 1;
model->setQuadratureRule(nGauss);
utl::profiler->stop("Model input");
// Write VTF-file with model geometry
int nBlock = 0;
if (!model->writeGlv(infile,n,format))
return 2;
// Write Dirichlet boundary conditions
if (!model->writeGlvBC(n,nBlock))
return 3;
Matrix eNorm;
Vector gNorm, displ, load;
model->initSystem(solver,1,1);
model->setAssociatedRHS(0,0);
for (int iStep = 1; iStep <= nStep; iStep++)
{
time.t += time.dt;
std::cout <<"\n\nStep "<< iStep <<", time = "<< time.t << std::endl;
model->setMode(SIM::STATIC);
if (!model->assembleSystem(time))
return 4;
else if (vizRHS)
model->extractLoadVec(load);
// Solve the linear system of equations
if (!model->solveSystem(displ,1))
return 5;
// Evaluate solution norms
model->setMode(SIM::RECOVERY);
if (!model->solutionNorms(Vectors(1,displ),eNorm,gNorm))
return 6;
std::cout <<"Energy norm |u^h| = a(u^h,u^h)^0.5 : "<< gNorm(1);
if (gNorm.size() > 1)
std::cout <<"\nExternal energy ((f,u^h)+(t,u^h)^0.5 : "<< gNorm(2);
if (gNorm.size() > 2)
std::cout <<"\nExact norm |u| = a(u,u)^0.5 : "<< gNorm(3);
if (gNorm.size() > 3)
std::cout <<"\nExact error a(e,e)^0.5, e=u-u^h : "<< gNorm(4)
<<"\nExact relative error (%) : "<< gNorm(4)/gNorm(3)*100.0;
std::cout << std::endl;
utl::profiler->start("Postprocessing");
// Write boundary tractions, if any
if (!model->writeGlvT(iStep,nBlock))
return 7;
// Write load vector to VTF-file
if (!model->writeGlvV(load,"Load vector",n,iStep,nBlock))
return 8;
// Write solution fields to VTF-file
if (!model->writeGlvS(displ,n,iStep,nBlock))
return 9;
// Write element norms, if feasable
if (n[0] == 2 && n[1] == 2 && n[2] == 2)
if (!model->writeGlvN(eNorm,iStep,nBlock))
return 10;
// Write step information to VTF-file
if (!model->writeGlvStep(iStep,time.t))
return 11;
utl::profiler->stop("Postprocessing");
}
model->closeGlv();
delete model;
return 0;
}