From cd5407b08b58f549f65575205d1964c8f506b87f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 25 Feb 2011 09:25:53 +0000 Subject: [PATCH] removed all math functions only for double precision by the more flexible counterpart, e.g. "dsqrt --> sqrt", "dsin --> sin". Should not cause any harm, as long as "implicit none" is used. Now it is possible to compile a single precision spectral solver/crystal plasticity by replacing mesh.f90 and prec.f90 with mesh_single.f90 and prec_single.f90. For the spectral method, just call "make precision=single" instead of "make". Use "make clean" evertime you switch precision --- code/constitutive_j2.f90 | 12 +- code/constitutive_nonlocal.f90 | 10 +- code/constitutive_phenopowerlaw.f90 | 2 +- code/lattice.f90 | 36 +- code/makefile | 58 +- code/makefile_single | 86 - code/math.f90 | 94 +- code/mesh.f90 | 2 +- code/mesh_single.f90 | 3360 +++++++++++++++++++++++ code/mpie_spectral_single.f90 | 2 +- processing/post/spectral_post_single.py | 368 +++ 11 files changed, 3859 insertions(+), 171 deletions(-) delete mode 100644 code/makefile_single create mode 100644 code/mesh_single.f90 create mode 100644 processing/post/spectral_post_single.py diff --git a/code/constitutive_j2.f90 b/code/constitutive_j2.f90 index 5e978c658..34b33a6ba 100644 --- a/code/constitutive_j2.f90 +++ b/code/constitutive_j2.f90 @@ -341,7 +341,7 @@ pure subroutine constitutive_j2_LpAndItsTangent(Lp, dLp_dTstar_99, Tstar_dev_v, ! convert Tstar to matrix and calculate euclidean norm Tstar_dev_33 = math_Mandel6to33(Tstar_dev_v) squarenorm_Tstar_dev = math_mul6x6(Tstar_dev_v,Tstar_dev_v) - norm_Tstar_dev = dsqrt(squarenorm_Tstar_dev) + norm_Tstar_dev = sqrt(squarenorm_Tstar_dev) ! Initialization of Lp and dLp_dTstar Lp = 0.0_pReal @@ -351,7 +351,7 @@ pure subroutine constitutive_j2_LpAndItsTangent(Lp, dLp_dTstar_99, Tstar_dev_v, if (norm_Tstar_dev > 0) then ! Calculation of gamma_dot - gamma_dot = constitutive_j2_gdot0(matID) * ( dsqrt(1.5_pReal) * norm_Tstar_dev & + gamma_dot = constitutive_j2_gdot0(matID) * ( sqrt(1.5_pReal) * norm_Tstar_dev & / &!--------------------------------------------------- (constitutive_j2_fTaylor(matID) * state(g,ip,el)%p(1)) ) **constitutive_j2_n(matID) @@ -414,10 +414,10 @@ pure function constitutive_j2_dotState(Tstar_v, Temperature, state, g, ip, el) Tstar_dev_v(1:3) = Tstar_v(1:3) - sum(Tstar_v(1:3))/3.0_pReal Tstar_dev_v(4:6) = Tstar_v(4:6) - norm_Tstar_dev = dsqrt(math_mul6x6(Tstar_dev_v,Tstar_dev_v)) + norm_Tstar_dev = sqrt(math_mul6x6(Tstar_dev_v,Tstar_dev_v)) ! gamma_dot - gamma_dot = constitutive_j2_gdot0(matID) * ( dsqrt(1.5_pReal) * norm_Tstar_dev & + gamma_dot = constitutive_j2_gdot0(matID) * ( sqrt(1.5_pReal) * norm_Tstar_dev & / &!--------------------------------------------------- (constitutive_j2_fTaylor(matID) * state(g,ip,el)%p(1)) ) ** constitutive_j2_n(matID) @@ -514,7 +514,7 @@ pure function constitutive_j2_postResults(Tstar_v, Temperature, dt, state, g, ip ! calculate deviatoric part of 2nd Piola-Kirchhoff stress and its norm Tstar_dev_v(1:3) = Tstar_v(1:3) - sum(Tstar_v(1:3))/3.0_pReal Tstar_dev_v(4:6) = Tstar_v(4:6) - norm_Tstar_dev = dsqrt(math_mul6x6(Tstar_dev_v,Tstar_dev_v)) + norm_Tstar_dev = sqrt(math_mul6x6(Tstar_dev_v,Tstar_dev_v)) c = 0_pInt constitutive_j2_postResults = 0.0_pReal @@ -526,7 +526,7 @@ pure function constitutive_j2_postResults(Tstar_v, Temperature, dt, state, g, ip c = c + 1 case ('strainrate') constitutive_j2_postResults(c+1) = & - constitutive_j2_gdot0(matID) * ( dsqrt(1.5_pReal) * norm_Tstar_dev & + constitutive_j2_gdot0(matID) * ( sqrt(1.5_pReal) * norm_Tstar_dev & / &!--------------------------------------------------- (constitutive_j2_fTaylor(matID) * state(g,ip,el)%p(1)) ) ** constitutive_j2_n(matID) c = c + 1 diff --git a/code/constitutive_nonlocal.f90 b/code/constitutive_nonlocal.f90 index e0c98ffda..46acc886d 100644 --- a/code/constitutive_nonlocal.f90 +++ b/code/constitutive_nonlocal.f90 @@ -636,7 +636,7 @@ do i = 1,maxNinstance !*** calculation of prefactor for activation enthalpy for dislocation glide constitutive_nonlocal_Qeff0(1:ns,i) = constitutive_nonlocal_burgersPerSlipSystem(1:ns,i) ** 3.0_pReal & - * dsqrt(0.5_pReal * constitutive_nonlocal_d0(i) ** 3.0_pReal & + * sqrt(0.5_pReal * constitutive_nonlocal_d0(i) ** 3.0_pReal & * constitutive_nonlocal_Gmod(i) * constitutive_nonlocal_tauObs(i)) enddo @@ -1122,10 +1122,10 @@ if (Temperature > 0.0_pReal) then tauRel = (abs(tau(s)) - tauThreshold(s)) / constitutive_nonlocal_tauObs(myInstance) if (tauRel > 0.0_pReal .and. tauRel < 1.0_pReal) then - wallFunc = 4.0_pReal * dsqrt(2.0_pReal) / 3.0_pReal * dsqrt(1.0_pReal - tauRel) / tauRel - boltzmannProbability = dexp(- constitutive_nonlocal_Qeff0(s,myInstance) * wallFunc / (kB * Temperature)) + wallFunc = 4.0_pReal * sqrt(2.0_pReal) / 3.0_pReal * sqrt(1.0_pReal - tauRel) / tauRel + boltzmannProbability = exp(- constitutive_nonlocal_Qeff0(s,myInstance) * wallFunc / (kB * Temperature)) timeRatio = boltzmannProbability * constitutive_nonlocal_fattack(myInstance) & - / (constitutive_nonlocal_vs(myInstance) * dsqrt(rhoForest(s))) + / (constitutive_nonlocal_vs(myInstance) * sqrt(rhoForest(s))) constitutive_nonlocal_v(s,:,g,ip,el) = sign(constitutive_nonlocal_vs(myInstance),tau(s)) * timeRatio / (1.0_pReal + timeRatio) @@ -1142,7 +1142,7 @@ if (Temperature > 0.0_pReal) then elseif (tauRel >= 1.0_pReal) then constitutive_nonlocal_v(s,1:4,g,ip,el) = sign(constitutive_nonlocal_vs(myInstance), tau(s)) & * constitutive_nonlocal_fattack(myInstance) & - / (constitutive_nonlocal_vs(myInstance) * dsqrt(rhoForest(s)) & + / (constitutive_nonlocal_vs(myInstance) * sqrt(rhoForest(s)) & + constitutive_nonlocal_fattack(myInstance)) endif enddo diff --git a/code/constitutive_phenopowerlaw.f90 b/code/constitutive_phenopowerlaw.f90 index a00f21757..13a12212d 100644 --- a/code/constitutive_phenopowerlaw.f90 +++ b/code/constitutive_phenopowerlaw.f90 @@ -763,7 +763,7 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el !-- add system-dependent part and calculate dot gammas - ssat_offset = constitutive_phenopowerlaw_spr(matID)*dsqrt(state(ipc,ip,el)%p(index_F)) + ssat_offset = constitutive_phenopowerlaw_spr(matID)*sqrt(state(ipc,ip,el)%p(index_F)) j = 0_pInt do f = 1,lattice_maxNslipFamily ! loop over all slip families index_myFamily = sum(lattice_NslipSystem(1:f-1,structID)) ! at which index starts my family diff --git a/code/lattice.f90 b/code/lattice.f90 index 4804dd201..9db7c60d7 100644 --- a/code/lattice.f90 +++ b/code/lattice.f90 @@ -776,13 +776,13 @@ function lattice_initializeStructure(struct,CoverA) if (lattice_fcc_Nstructure == 1_pInt) then ! me is first fcc structure processMe = .true. do i = 1,myNslip ! calculate slip system vectors - sd(:,i) = lattice_fcc_systemSlip(1:3,i)/dsqrt(math_mul3x3(lattice_fcc_systemSlip(1:3,i),lattice_fcc_systemSlip(1:3,i))) - sn(:,i) = lattice_fcc_systemSlip(4:6,i)/dsqrt(math_mul3x3(lattice_fcc_systemSlip(4:6,i),lattice_fcc_systemSlip(4:6,i))) + sd(:,i) = lattice_fcc_systemSlip(1:3,i)/sqrt(math_mul3x3(lattice_fcc_systemSlip(1:3,i),lattice_fcc_systemSlip(1:3,i))) + sn(:,i) = lattice_fcc_systemSlip(4:6,i)/sqrt(math_mul3x3(lattice_fcc_systemSlip(4:6,i),lattice_fcc_systemSlip(4:6,i))) st(:,i) = math_vectorproduct(sd(:,i),sn(:,i)) enddo do i = 1,myNtwin ! calculate twin system vectors and (assign) shears - td(:,i) = lattice_fcc_systemTwin(1:3,i)/dsqrt(math_mul3x3(lattice_fcc_systemTwin(1:3,i),lattice_fcc_systemTwin(1:3,i))) - tn(:,i) = lattice_fcc_systemTwin(4:6,i)/dsqrt(math_mul3x3(lattice_fcc_systemTwin(4:6,i),lattice_fcc_systemTwin(4:6,i))) + td(:,i) = lattice_fcc_systemTwin(1:3,i)/sqrt(math_mul3x3(lattice_fcc_systemTwin(1:3,i),lattice_fcc_systemTwin(1:3,i))) + tn(:,i) = lattice_fcc_systemTwin(4:6,i)/sqrt(math_mul3x3(lattice_fcc_systemTwin(4:6,i),lattice_fcc_systemTwin(4:6,i))) tt(:,i) = math_vectorproduct(td(:,i),tn(:,i)) ts(i) = lattice_fcc_shearTwin(i) enddo @@ -802,13 +802,13 @@ function lattice_initializeStructure(struct,CoverA) if (lattice_bcc_Nstructure == 1_pInt) then ! me is first bcc structure processMe = .true. do i = 1,myNslip ! calculate slip system vectors - sd(:,i) = lattice_bcc_systemSlip(1:3,i)/dsqrt(math_mul3x3(lattice_bcc_systemSlip(1:3,i),lattice_bcc_systemSlip(1:3,i))) - sn(:,i) = lattice_bcc_systemSlip(4:6,i)/dsqrt(math_mul3x3(lattice_bcc_systemSlip(4:6,i),lattice_bcc_systemSlip(4:6,i))) + sd(:,i) = lattice_bcc_systemSlip(1:3,i)/sqrt(math_mul3x3(lattice_bcc_systemSlip(1:3,i),lattice_bcc_systemSlip(1:3,i))) + sn(:,i) = lattice_bcc_systemSlip(4:6,i)/sqrt(math_mul3x3(lattice_bcc_systemSlip(4:6,i),lattice_bcc_systemSlip(4:6,i))) st(:,i) = math_vectorproduct(sd(:,i),sn(:,i)) enddo do i = 1,myNtwin ! calculate twin system vectors and (assign) shears - td(:,i) = lattice_bcc_systemTwin(1:3,i)/dsqrt(math_mul3x3(lattice_bcc_systemTwin(1:3,i),lattice_bcc_systemTwin(1:3,i))) - tn(:,i) = lattice_bcc_systemTwin(4:6,i)/dsqrt(math_mul3x3(lattice_bcc_systemTwin(4:6,i),lattice_bcc_systemTwin(4:6,i))) + td(:,i) = lattice_bcc_systemTwin(1:3,i)/sqrt(math_mul3x3(lattice_bcc_systemTwin(1:3,i),lattice_bcc_systemTwin(1:3,i))) + tn(:,i) = lattice_bcc_systemTwin(4:6,i)/sqrt(math_mul3x3(lattice_bcc_systemTwin(4:6,i),lattice_bcc_systemTwin(4:6,i))) tt(:,i) = math_vectorproduct(td(:,i),tn(:,i)) ts(i) = lattice_bcc_shearTwin(i) enddo @@ -830,37 +830,37 @@ function lattice_initializeStructure(struct,CoverA) ! converting from 4 axes coordinate system (a1=a2=a3=c) to ortho-hexgonal system (a, b, c) do i = 1,myNslip hex_d(1) = lattice_hex_systemSlip(1,i)*1.5_pReal ! direction [uvtw]->[3u/2 (u+2v)*sqrt(3)/2 w*(c/a)] - hex_d(2) = (lattice_hex_systemSlip(1,i)+2.0_pReal*lattice_hex_systemSlip(2,i))*(0.5_pReal*dsqrt(3.0_pReal)) + hex_d(2) = (lattice_hex_systemSlip(1,i)+2.0_pReal*lattice_hex_systemSlip(2,i))*(0.5_pReal*sqrt(3.0_pReal)) hex_d(3) = lattice_hex_systemSlip(4,i)*CoverA hex_n(1) = lattice_hex_systemSlip(5,i) ! plane (hkil)->(h (h+2k)/sqrt(3) l/(c/a)) - hex_n(2) = (lattice_hex_systemSlip(5,i)+2.0_pReal*lattice_hex_systemSlip(6,i))/dsqrt(3.0_pReal) + hex_n(2) = (lattice_hex_systemSlip(5,i)+2.0_pReal*lattice_hex_systemSlip(6,i))/sqrt(3.0_pReal) hex_n(3) = lattice_hex_systemSlip(8,i)/CoverA - sd(:,i) = hex_d/dsqrt(math_mul3x3(hex_d,hex_d)) - sn(:,i) = hex_n/dsqrt(math_mul3x3(hex_n,hex_n)) + sd(:,i) = hex_d/sqrt(math_mul3x3(hex_d,hex_d)) + sn(:,i) = hex_n/sqrt(math_mul3x3(hex_n,hex_n)) st(:,i) = math_vectorproduct(sd(:,i),sn(:,i)) enddo do i = 1,myNtwin hex_d(1) = lattice_hex_systemTwin(1,i)*1.5_pReal - hex_d(2) = (lattice_hex_systemTwin(1,i)+2.0_pReal*lattice_hex_systemTwin(2,i))*(0.5_pReal*dsqrt(3.0_pReal)) + hex_d(2) = (lattice_hex_systemTwin(1,i)+2.0_pReal*lattice_hex_systemTwin(2,i))*(0.5_pReal*sqrt(3.0_pReal)) hex_d(3) = lattice_hex_systemTwin(4,i)*CoverA hex_n(1) = lattice_hex_systemTwin(5,i) - hex_n(2) = (lattice_hex_systemTwin(5,i)+2.0_pReal*lattice_hex_systemTwin(6,i))/dsqrt(3.0_pReal) + hex_n(2) = (lattice_hex_systemTwin(5,i)+2.0_pReal*lattice_hex_systemTwin(6,i))/sqrt(3.0_pReal) hex_n(3) = lattice_hex_systemTwin(8,i)/CoverA - td(:,i) = hex_d/dsqrt(math_mul3x3(hex_d,hex_d)) - tn(:,i) = hex_n/dsqrt(math_mul3x3(hex_n,hex_n)) + td(:,i) = hex_d/sqrt(math_mul3x3(hex_d,hex_d)) + tn(:,i) = hex_n/sqrt(math_mul3x3(hex_n,hex_n)) tt(:,i) = math_vectorproduct(td(:,i),tn(:,i)) select case(lattice_hex_shearTwin(i)) ! from Christian & Mahajan 1995 p.29 case (1) ! {10.2}<-10.1> - ts(i) = (3.0_pReal-CoverA*CoverA)/dsqrt(3.0_pReal)/CoverA + ts(i) = (3.0_pReal-CoverA*CoverA)/sqrt(3.0_pReal)/CoverA case (2) ! {11.2}<11.-3> ts(i) = 2.0_pReal*(CoverA*CoverA-2.0_pReal)/3.0_pReal/CoverA case (3) ! {11.1}<-1-1.6> ts(i) = 1.0_pReal/CoverA case (4) ! {10.1}<10.-2> - ts(i) = (4.0_pReal*CoverA*CoverA-9.0_pReal)/4.0_pReal/dsqrt(3.0_pReal)/CoverA + ts(i) = (4.0_pReal*CoverA*CoverA-9.0_pReal)/4.0_pReal/sqrt(3.0_pReal)/CoverA end select enddo diff --git a/code/makefile b/code/makefile index f8ce20266..01bc539bb 100644 --- a/code/makefile +++ b/code/makefile @@ -1,17 +1,32 @@ -# Makefile to compile the Material subroutine for BVP solution using spectral method. +# Makefile to compile the Material subroutine for BVP solution using spectral method # +# use switch on make to determine precision, e.g make precision=single +# default is precision=double +# be sure to remove all librarys with different precision (make clean) +# # Uses openmp to parallelise the material subroutines (set number of cores with "export MPIE_NUM_THREADS=n" to n) # Uses linux threads to parallelise fftw3 (should also be possible with openmp) # Besides of the f90 files written at MPIE, the two library files of fftw3 "libfftw3_threads.a" "libfftw3.a" are also needed -# Install fftw3 (v3.2.2 is tested) with "./configure --enable-threads" and "make", "make install" is not needed +# Install fftw3 (v3.2.2 is tested) with "./configure --enable-threads --enable-float" and "make", "make install" is not needed # as long as the two library files are copied to the source code directory. -cpspectral.exe: mpie_spectral.o CPFEM.a - ifort -openmp -o cpspectral.exe mpie_spectral.o CPFEM.a libfftw3_threads.a libfftw3.a constitutive.a advanced.a basics.a -lpthread +precision :=double +ifeq ($(precision),single) +cpspectral_single.exe: mpie_spectral_single.o CPFEM.a + ifort -openmp -o cpspectral_single.exe mpie_spectral_single.o CPFEM.a libfftw3f_threads.a libfftw3f.a constitutive.a advanced.a basics.a -lpthread + rm *.mod + +mpie_spectral_single.o: mpie_spectral_single.f90 CPFEM.o + ifort -openmp -c -O3 -heap-arrays 500000000 mpie_spectral_single.f90 + +else +cpspectral_double.exe: mpie_spectral.o CPFEM.a + ifort -openmp -o cpspectral_double.exe mpie_spectral.o CPFEM.a libfftw3_threads.a libfftw3.a constitutive.a advanced.a basics.a -lpthread rm *.mod mpie_spectral.o: mpie_spectral.f90 CPFEM.o ifort -openmp -c -O3 -heap-arrays 500000000 mpie_spectral.f90 +endif CPFEM.a: CPFEM.o ar rc CPFEM.a homogenization.o homogenization_RGC.o homogenization_isostrain.o crystallite.o CPFEM.o constitutive.o @@ -51,23 +66,38 @@ constitutive_phenopowerlaw.o: constitutive_phenopowerlaw.f90 basics.a advanced.a ifort -openmp -c -O3 -heap-arrays 500000000 constitutive_phenopowerlaw.f90 - +ifeq ($(precision),single) +advanced.a: lattice.o + ar rc advanced.a FEsolving.o mesh_single.o material.o lattice.o +else advanced.a: lattice.o ar rc advanced.a FEsolving.o mesh.o material.o lattice.o +endif lattice.o: lattice.f90 material.o ifort -openmp -c -O3 -heap-arrays 500000000 lattice.f90 +ifeq ($(precision),single) +material.o: material.f90 mesh_single.o + ifort -openmp -c -O3 -heap-arrays 500000000 material.f90 +mesh_single.o: mesh_single.f90 FEsolving.o + ifort -openmp -c -O3 -heap-arrays 500000000 mesh_single.f90 +else material.o: material.f90 mesh.o ifort -openmp -c -O3 -heap-arrays 500000000 material.f90 mesh.o: mesh.f90 FEsolving.o ifort -openmp -c -O3 -heap-arrays 500000000 mesh.f90 +endif FEsolving.o: FEsolving.f90 basics.a ifort -openmp -c -O3 -heap-arrays 500000000 FEsolving.f90 - +ifeq ($(precision),single) +basics.a: debug.o math.o + ar rc basics.a debug.o math.o numerics.o IO.o mpie_spectral_interface.o prec_single.o +else basics.a: debug.o math.o ar rc basics.a debug.o math.o numerics.o IO.o mpie_spectral_interface.o prec.o +endif debug.o: debug.f90 numerics.o ifort -openmp -c -O3 debug.f90 @@ -79,7 +109,23 @@ numerics.o: numerics.f90 IO.o ifort -openmp -c -O3 numerics.f90 IO.o: IO.f90 mpie_spectral_interface.o ifort -openmp -c -O3 IO.f90 +ifeq ($(precision),single) +mpie_spectral_interface.o: mpie_spectral_interface.f90 prec_single.o + ifort -openmp -c -O3 mpie_spectral_interface.f90 +prec_single.o: prec_single.f90 + ifort -openmp -c -O3 prec_single.f90 +else mpie_spectral_interface.o: mpie_spectral_interface.f90 prec.o ifort -openmp -c -O3 mpie_spectral_interface.f90 prec.o: prec.f90 ifort -openmp -c -O3 prec.f90 +endif + + +clean: + rm -rf *.o + rm -rf *.exe + rm -rf basics.a + rm -rf advanced.a + rm -rf constitutive.a + rm -rf CPFEM.a diff --git a/code/makefile_single b/code/makefile_single deleted file mode 100644 index ea2e1c071..000000000 --- a/code/makefile_single +++ /dev/null @@ -1,86 +0,0 @@ -# Makefile to compile the Material subroutine for BVP solution using spectral method with single precision -# -# BE SURE THAT ALL DOUBLE FILES ARE REMOVE FROM THE FOLDER -# Uses openmp to parallelise the material subroutines (set number of cores with "export MPIE_NUM_THREADS=n" to n) -# Uses linux threads to parallelise fftw3 (should also be possible with openmp) -# Besides of the f90 files written at MPIE, the two library files of fftw3 "libfftw3_threads.a" "libfftw3.a" are also needed -# Install fftw3 (v3.2.2 is tested) with "./configure --enable-threads --enable-float" and "make", "make install" is not needed -# as long as the two library files are copied to the source code directory. - -cpspectral_single.exe: mpie_spectral.o CPFEM.a - ifort -openmp -o cpspectral_single.exe mpie_spectral.o CPFEM.a libfftw3f_threads.a libfftw3f.a constitutive.a advanced.a basics.a -lpthread - rm *.mod - -mpie_spectral.o: mpie_spectral.f90 CPFEM.o - ifort -openmp -c -O3 -heap-arrays 500000000 mpie_spectral.f90 - -CPFEM.a: CPFEM.o - ar rc CPFEM.a homogenization.o homogenization_RGC.o homogenization_isostrain.o crystallite.o CPFEM.o constitutive.o - -CPFEM.o: CPFEM.f90 homogenization.o - ifort -openmp -c -O3 -heap-arrays 500000000 CPFEM.f90 -homogenization.o: homogenization.f90 homogenization_isostrain.o homogenization_RGC.o crystallite.o - ifort -openmp -c -O3 -heap-arrays 500000000 homogenization.f90 -homogenization_RGC.o: homogenization_RGC.f90 constitutive.a - ifort -openmp -c -O3 -heap-arrays 500000000 homogenization_RGC.f90 -homogenization_isostrain.o: homogenization_isostrain.f90 basics.a advanced.a - ifort -openmp -c -O3 -heap-arrays 500000000 homogenization_isostrain.f90 -crystallite.o: crystallite.f90 constitutive.a - ifort -openmp -c -O3 -heap-arrays 500000000 crystallite.f90 - - - -constitutive.a: constitutive.o - ar rc constitutive.a constitutive.o constitutive_titanmod.o constitutive_nonlocal.o constitutive_dislotwin.o constitutive_j2.o constitutive_phenopowerlaw.o basics.a advanced.a - -constitutive.o: constitutive.f90 constitutive_titanmod.o constitutive_nonlocal.o constitutive_dislotwin.o constitutive_j2.o constitutive_phenopowerlaw.o - ifort -openmp -c -O3 -heap-arrays 500000000 constitutive.f90 - -constitutive_titanmod.o: constitutive_titanmod.f90 basics.a advanced.a - ifort -openmp -c -O3 -heap-arrays 500000000 constitutive_titanmod.f90 - -constitutive_nonlocal.o: constitutive_nonlocal.f90 basics.a advanced.a - ifort -openmp -c -O3 -heap-arrays 500000000 constitutive_nonlocal.f90 - -constitutive_dislotwin.o: constitutive_dislotwin.f90 basics.a advanced.a - ifort -openmp -c -O3 -heap-arrays 500000000 constitutive_dislotwin.f90 - -constitutive_j2.o: constitutive_j2.f90 basics.a advanced.a - ifort -openmp -c -O3 -heap-arrays 500000000 constitutive_j2.f90 - -constitutive_phenopowerlaw.o: constitutive_phenopowerlaw.f90 basics.a advanced.a - ifort -openmp -c -O3 -heap-arrays 500000000 constitutive_phenopowerlaw.f90 - - - -advanced.a: lattice.o - ar rc advanced.a FEsolving.o mesh.o material.o lattice.o - -lattice.o: lattice.f90 material.o - ifort -openmp -c -O3 -heap-arrays 500000000 lattice.f90 -material.o: material.f90 mesh.o - ifort -openmp -c -O3 -heap-arrays 500000000 material.f90 -mesh.o: mesh.f90 FEsolving.o - ifort -openmp -c -O3 -heap-arrays 500000000 mesh.f90 -FEsolving.o: FEsolving.f90 basics.a - ifort -openmp -c -O3 -heap-arrays 500000000 FEsolving.f90 - - - -basics.a: debug.o math_single.o - ar rc basics.a debug.o math_single.o numerics.o IO.o mpie_spectral_interface.o prec_single.o - -debug.o: debug.f90 numerics.o - ifort -openmp -c -O3 debug.f90 - -math_single.o: math_single.f90 numerics.o - ifort -openmp -c -O3 math_single.f90 - -numerics.o: numerics.f90 IO.o - ifort -openmp -c -O3 numerics.f90 -IO.o: IO.f90 mpie_spectral_interface.o - ifort -openmp -c -O3 IO.f90 -mpie_spectral_interface.o: mpie_spectral_interface.f90 prec_single.o - ifort -openmp -c -O3 mpie_spectral_interface.f90 -prec_single.o: prec_single.f90 - ifort -openmp -c -O3 prec_single.f90 diff --git a/code/math.f90 b/code/math.f90 index ba9fbb297..401f28b7d 100644 --- a/code/math.f90 +++ b/code/math.f90 @@ -604,10 +604,10 @@ real(pReal), dimension(4,36), parameter :: math_symOperations = & real(pReal), dimension(3) :: rnd call halton(3,rnd) - math_qRnd(1) = dcos(2.0_pReal*pi*rnd(1))*dsqrt(rnd(3)) - math_qRnd(2) = dsin(2.0_pReal*pi*rnd(2))*dsqrt(1.0_pReal-rnd(3)) - math_qRnd(3) = dcos(2.0_pReal*pi*rnd(2))*dsqrt(1.0_pReal-rnd(3)) - math_qRnd(4) = dsin(2.0_pReal*pi*rnd(1))*dsqrt(rnd(3)) + math_qRnd(1) = cos(2.0_pReal*pi*rnd(1))*sqrt(rnd(3)) + math_qRnd(2) = sin(2.0_pReal*pi*rnd(2))*sqrt(1.0_pReal-rnd(3)) + math_qRnd(3) = cos(2.0_pReal*pi*rnd(2))*sqrt(1.0_pReal-rnd(3)) + math_qRnd(4) = sin(2.0_pReal*pi*rnd(1))*sqrt(rnd(3)) endfunction @@ -675,7 +675,7 @@ real(pReal), dimension(4,36), parameter :: math_symOperations = & real(pReal), dimension(4), intent(in) :: Q real(pReal) math_qNorm - math_qNorm = dsqrt(max(0.0_pReal, Q(1)*Q(1) + Q(2)*Q(2) + Q(3)*Q(3) + Q(4)*Q(4))) + math_qNorm = sqrt(max(0.0_pReal, Q(1)*Q(1) + Q(2)*Q(2) + Q(3)*Q(3) + Q(4)*Q(4))) endfunction @@ -1127,7 +1127,7 @@ pure function math_transpose3x3(A) real(pReal), dimension(3), intent(in) :: v real(pReal) math_norm3 - math_norm3 = dsqrt(v(1)*v(1) + v(2)*v(2) + v(3)*v(3)) + math_norm3 = sqrt(v(1)*v(1) + v(2)*v(2) + v(3)*v(3)) return endfunction @@ -1333,7 +1333,7 @@ pure function math_transpose3x3(A) if(val > 1.0_pReal) val = 1.0_pReal if(val < -1.0_pReal) val = -1.0_pReal - math_RtoEuler(2) = dacos(val) + math_RtoEuler(2) = acos(val) if(math_RtoEuler(2) < 1.0e-8_pReal) then ! calculate phi2 @@ -1343,7 +1343,7 @@ pure function math_transpose3x3(A) if(val > 1.0_pReal) val = 1.0_pReal if(val < -1.0_pReal) val = -1.0_pReal - math_RtoEuler(1) = dacos(val) + math_RtoEuler(1) = acos(val) if(R(2,1) > 0.0_pReal) math_RtoEuler(1) = 2.0_pReal*pi-math_RtoEuler(1) else ! calculate phi2 @@ -1351,14 +1351,14 @@ pure function math_transpose3x3(A) if(val > 1.0_pReal) val = 1.0_pReal if(val < -1.0_pReal) val = -1.0_pReal - math_RtoEuler(3) = dacos(val) + math_RtoEuler(3) = acos(val) if(R(1,3) < 0.0) math_RtoEuler(3) = 2.0_pReal*pi-math_RtoEuler(3) ! calculate phi1 val=-R(3,2)/sin(math_RtoEuler(2)) if(val > 1.0_pReal) val = 1.0_pReal if(val < -1.0_pReal) val = -1.0_pReal - math_RtoEuler(1) = dacos(val) + math_RtoEuler(1) = acos(val) if(R(3,1) < 0.0) math_RtoEuler(1) = 2.0_pReal*pi-math_RtoEuler(1) end if return @@ -1382,10 +1382,10 @@ pure function math_transpose3x3(A) math_RtoQuaternion = 0.0_pReal - absQ(1) = 0.5_pReal * dsqrt(1.0_pReal+R(1,1)+R(2,2)+R(3,3)) - absQ(2) = 0.5_pReal * dsqrt(1.0_pReal+R(1,1)-R(2,2)-R(3,3)) - absQ(3) = 0.5_pReal * dsqrt(1.0_pReal-R(1,1)+R(2,2)-R(3,3)) - absQ(4) = 0.5_pReal * dsqrt(1.0_pReal-R(1,1)-R(2,2)+R(3,3)) + absQ(1) = 0.5_pReal * sqrt(1.0_pReal+R(1,1)+R(2,2)+R(3,3)) + absQ(2) = 0.5_pReal * sqrt(1.0_pReal+R(1,1)-R(2,2)-R(3,3)) + absQ(3) = 0.5_pReal * sqrt(1.0_pReal-R(1,1)+R(2,2)-R(3,3)) + absQ(4) = 0.5_pReal * sqrt(1.0_pReal-R(1,1)-R(2,2)+R(3,3)) largest = maxloc(absQ) select case(largest(1)) @@ -1434,12 +1434,12 @@ pure function math_transpose3x3(A) real(pReal), dimension(3,3) :: math_EulerToR real(pReal) c1, c, c2, s1, s, s2 - C1 = dcos(Euler(1)) - C = dcos(Euler(2)) - C2 = dcos(Euler(3)) - S1 = dsin(Euler(1)) - S = dsin(Euler(2)) - S2 = dsin(Euler(3)) + C1 = cos(Euler(1)) + C = cos(Euler(2)) + C2 = cos(Euler(3)) + S1 = sin(Euler(1)) + S = sin(Euler(2)) + S2 = sin(Euler(3)) math_EulerToR(1,1)=C1*C2-S1*S2*C math_EulerToR(1,2)=S1*C2+C1*S2*C math_EulerToR(1,3)=S2*S @@ -1469,13 +1469,13 @@ pure function math_transpose3x3(A) halfangles = 0.5_pReal * eulerangles - c = dcos(halfangles(2)) - s = dsin(halfangles(2)) + c = cos(halfangles(2)) + s = sin(halfangles(2)) - math_EulerToQuaternion(1) = dcos(halfangles(1)+halfangles(3)) * c - math_EulerToQuaternion(2) = dcos(halfangles(1)-halfangles(3)) * s - math_EulerToQuaternion(3) = dsin(halfangles(1)-halfangles(3)) * s - math_EulerToQuaternion(4) = dsin(halfangles(1)+halfangles(3)) * c + math_EulerToQuaternion(1) = cos(halfangles(1)+halfangles(3)) * c + math_EulerToQuaternion(2) = cos(halfangles(1)-halfangles(3)) * s + math_EulerToQuaternion(3) = sin(halfangles(1)-halfangles(3)) * s + math_EulerToQuaternion(4) = sin(halfangles(1)+halfangles(3)) * c endfunction @@ -1495,12 +1495,12 @@ pure function math_transpose3x3(A) real(pReal) norm,s,c,c1 integer(pInt) i - norm = dsqrt(math_mul3x3(axis,axis)) + norm = sqrt(math_mul3x3(axis,axis)) if (norm > 1.0e-8_pReal) then ! non-zero rotation forall (i=1:3) axisNrm(i) = axis(i)/norm ! normalize axis to be sure - s = dsin(omega) - c = dcos(omega) + s = sin(omega) + c = cos(omega) c1 = 1.0_pReal - c ! formula for active rotation taken from http://mathworld.wolfram.com/RodriguesRotationFormula.html @@ -1541,12 +1541,12 @@ pure function math_transpose3x3(A) real(pReal) s,c,norm integer(pInt) i - norm = dsqrt(math_mul3x3(axis,axis)) + norm = sqrt(math_mul3x3(axis,axis)) if (norm > 1.0e-8_pReal) then ! non-zero rotation forall (i=1:3) axisNrm(i) = axis(i)/norm ! normalize axis to be sure ! formula taken from http://en.wikipedia.org/wiki/Rotation_representation_%28mathematics%29#Rodrigues_parameters - s = dsin(omega/2.0_pReal) - c = dcos(omega/2.0_pReal) + s = sin(omega/2.0_pReal) + c = cos(omega/2.0_pReal) math_AxisAngleToQuaternion(1) = c math_AxisAngleToQuaternion(2:4) = s * axisNrm(1:3) else @@ -1597,17 +1597,17 @@ pure function math_transpose3x3(A) real(pReal), dimension(4), intent(in) :: Q real(pReal), dimension(3) :: math_QuaternionToEuler - math_QuaternionToEuler(2) = dacos(1.0_pReal-2.0_pReal*(Q(2)*Q(2)+Q(3)*Q(3))) + math_QuaternionToEuler(2) = acos(1.0_pReal-2.0_pReal*(Q(2)*Q(2)+Q(3)*Q(3))) - if (dabs(math_QuaternionToEuler(2)) < 1.0e-3_pReal) then - math_QuaternionToEuler(1) = 2.0_pReal*dacos(Q(1)) + if (abs(math_QuaternionToEuler(2)) < 1.0e-3_pReal) then + math_QuaternionToEuler(1) = 2.0_pReal*acos(Q(1)) math_QuaternionToEuler(3) = 0.0_pReal else - math_QuaternionToEuler(1) = datan2(Q(1)*Q(3)+Q(2)*Q(4), Q(1)*Q(2)-Q(3)*Q(4)) + math_QuaternionToEuler(1) = atan2(Q(1)*Q(3)+Q(2)*Q(4), Q(1)*Q(2)-Q(3)*Q(4)) if (math_QuaternionToEuler(1) < 0.0_pReal) & math_QuaternionToEuler(1) = math_QuaternionToEuler(1) + 2.0_pReal * pi - math_QuaternionToEuler(3) = datan2(-Q(1)*Q(3)+Q(2)*Q(4), Q(1)*Q(2)+Q(3)*Q(4)) + math_QuaternionToEuler(3) = atan2(-Q(1)*Q(3)+Q(2)*Q(4), Q(1)*Q(2)+Q(3)*Q(4)) if (math_QuaternionToEuler(3) < 0.0_pReal) & math_QuaternionToEuler(3) = math_QuaternionToEuler(3) + 2.0_pReal * pi endif @@ -1632,8 +1632,8 @@ pure function math_transpose3x3(A) real(pReal) halfAngle, sinHalfAngle real(pReal), dimension(4) :: math_QuaternionToAxisAngle - halfAngle = dacos(max(-1.0_pReal, min(1.0_pReal, Q(1)))) ! limit to [-1,1] --> 0 to 180 deg - sinHalfAngle = dsin(halfAngle) + halfAngle = acos(max(-1.0_pReal, min(1.0_pReal, Q(1)))) ! limit to [-1,1] --> 0 to 180 deg + sinHalfAngle = sin(halfAngle) if (sinHalfAngle <= 1.0e-4_pReal) then ! very small rotation angle? math_QuaternionToAxisAngle = 0.0_pReal @@ -1715,7 +1715,7 @@ pure function math_QuaternionInSST(Q, symmetryType) Rodrig(2) > Rodrig(3) .and. & Rodrig(3) > 0.0_pReal case (2) - math_QuaternionInSST = Rodrig(1) > dsqrt(3.0_pReal)*Rodrig(2) .and. & + math_QuaternionInSST = Rodrig(1) > sqrt(3.0_pReal)*Rodrig(2) .and. & Rodrig(2) > 0.0_pReal .and. & Rodrig(3) > 0.0_pReal case default @@ -1868,7 +1868,7 @@ endif fiberInS(3)=cos(beta(1)) ! ---# rotation matrix from sample to crystal system #--- - angle = -dacos(dot_product(fiberInC,fiberInS)) + angle = -acos(dot_product(fiberInC,fiberInS)) if(angle /= 0.0_pReal) then ! rotation axis between sample and crystal system (cross product) forall(i=1:3) axis(i) = fiberInC(rotMap(1,i))*fiberInS(rotMap(2,i))-fiberInC(rotMap(2,i))*fiberInS(rotMap(1,i)) @@ -2018,7 +2018,7 @@ endfunction ce = math_mul33x33(transpose(FE),FE) CALL math_spectral1(CE,EW1,EW2,EW3,EB1,EB2,EB3) - U=DSQRT(EW1)*EB1+DSQRT(EW2)*EB2+DSQRT(EW3)*EB3 + U=sqrt(EW1)*EB1+sqrt(EW2)*EB2+sqrt(EW3)*EB3 call math_invert3x3(U,UI,det,error) if (.not. error) R = math_mul33x33(FE,UI) @@ -2059,14 +2059,14 @@ endfunction EB2(2,2)=1.0_pReal EB3(3,3)=1.0_pReal ELSE - RHO=DSQRT(-3.0_pReal*P**3.0_pReal)/9.0_pReal + RHO=sqrt(-3.0_pReal*P**3.0_pReal)/9.0_pReal arg=-Q/RHO/2.0_pReal if(arg.GT.1) arg=1 if(arg.LT.-1) arg=-1 - PHI=DACOS(arg) - Y1=2*RHO**(1.0_pReal/3.0_pReal)*DCOS(PHI/3.0_pReal) - Y2=2*RHO**(1.0_pReal/3.0_pReal)*DCOS(PHI/3.0_pReal+2.0_pReal/3.0_pReal*PI) - Y3=2*RHO**(1.0_pReal/3.0_pReal)*DCOS(PHI/3.0_pReal+4.0_pReal/3.0_pReal*PI) + PHI=acos(arg) + Y1=2*RHO**(1.0_pReal/3.0_pReal)*cos(PHI/3.0_pReal) + Y2=2*RHO**(1.0_pReal/3.0_pReal)*cos(PHI/3.0_pReal+2.0_pReal/3.0_pReal*PI) + Y3=2*RHO**(1.0_pReal/3.0_pReal)*cos(PHI/3.0_pReal+4.0_pReal/3.0_pReal*PI) EW1=Y1-R/3.0_pReal EW2=Y2-R/3.0_pReal EW3=Y3-R/3.0_pReal diff --git a/code/mesh.f90 b/code/mesh.f90 index a72b35572..97f23716e 100644 --- a/code/mesh.f90 +++ b/code/mesh.f90 @@ -3127,7 +3127,7 @@ endsubroutine forall (n = 1:FE_NipFaceNodes, j = 1:Ntriangles) ! start at each interface node and build valid triangles to cover interface normal(:,j,n) = math_vectorproduct(nPos(:,1+mod(n+j-1,FE_NipFaceNodes)) - nPos(:,n), & ! calc their normal vectors nPos(:,1+mod(n+j-0,FE_NipFaceNodes)) - nPos(:,n)) - area(j,n) = dsqrt(sum(normal(:,j,n)*normal(:,j,n))) ! and area + area(j,n) = sqrt(sum(normal(:,j,n)*normal(:,j,n))) ! and area end forall forall (n = 1:FE_NipFaceNodes, j = 1:Ntriangles, area(j,n) > 0.0_pReal) & normal(:,j,n) = normal(:,j,n) / area(j,n) ! make unit normal diff --git a/code/mesh_single.f90 b/code/mesh_single.f90 new file mode 100644 index 000000000..ce4d7311f --- /dev/null +++ b/code/mesh_single.f90 @@ -0,0 +1,3360 @@ +!* $Id: mesh.f90 765 2011-02-16 16:23:08Z MPIE\c.kords $ +!############################################################## + MODULE mesh +!############################################################## + + use prec, only: pReal,pInt + implicit none + +! --------------------------- +! _Nelems : total number of elements in mesh +! _NcpElems : total number of CP elements in mesh +! _Nnodes : total number of nodes in mesh +! _maxNnodes : max number of nodes in any CP element +! _maxNips : max number of IPs in any CP element +! _maxNipNeighbors : max number of IP neighbors in any CP element +! _maxNsharedElems : max number of CP elements sharing a node +! +! _element : FEid, type(internal representation), material, texture, node indices +! _node : x,y,z coordinates (initially!) +! _sharedElem : entryCount and list of elements containing node +! +! _mapFEtoCPelem : [sorted FEid, corresponding CPid] +! _mapFEtoCPnode : [sorted FEid, corresponding CPid] +! +! MISSING: these definitions should actually reside in the +! FE-solver specific part (different for MARC/ABAQUS)..! +! Hence, I suggest to prefix with "FE_" +! +! _Nnodes : # nodes in a specific type of element (how we use it) +! _NoriginalNodes : # nodes in a specific type of element (how it is originally defined by marc) +! _Nips : # IPs in a specific type of element +! _NipNeighbors : # IP neighbors in a specific type of element +! _ipNeighbor : +x,-x,+y,-y,+z,-z list of intra-element IPs and +! (negative) neighbor faces per own IP in a specific type of element +! _NfaceNodes : # nodes per face in a specific type of element + +! _nodeOnFace : list of node indices on each face of a specific type of element +! _maxNnodesAtIP : max number of (equivalent) nodes attached to an IP +! _nodesAtIP : map IP index to two node indices in a specific type of element +! _ipNeighborhood : 6 or less neighboring IPs as [element_num, IP_index] +! _NsubNodes : # subnodes required to fully define all IP volumes + +! order is +x,-x,+y,-y,+z,-z but meaning strongly depends on Elemtype +! --------------------------- + integer(pInt) mesh_Nelems, mesh_NcpElems, mesh_NelemSets, mesh_maxNelemInSet + integer(pInt) mesh_Nmaterials + integer(pInt) mesh_Nnodes, mesh_maxNnodes, mesh_maxNips, mesh_maxNipNeighbors, mesh_maxNsharedElems, mesh_maxNsubNodes + integer(pInt), dimension(2) :: mesh_maxValStateVar = 0_pInt + character(len=64), dimension(:), allocatable :: mesh_nameElemSet, & ! names of elementSet + mesh_nameMaterial, & ! names of material in solid section + mesh_mapMaterial ! name of elementSet for material + integer(pInt), dimension(:,:), allocatable :: mesh_mapElemSet ! list of elements in elementSet + integer(pInt), dimension(:,:), allocatable, target :: mesh_mapFEtoCPelem, mesh_mapFEtoCPnode + integer(pInt), dimension(:,:), allocatable :: mesh_element, & + mesh_sharedElem, & + mesh_nodeTwins ! node twins are surface nodes that lie exactly on opposite sides of the mesh (surfaces nodes with equal coordinate values in two dimensions) + integer(pInt), dimension(:,:,:,:), allocatable :: mesh_ipNeighborhood + + real(pReal), dimension(:,:,:), allocatable :: mesh_subNodeCoord ! coordinates of subnodes per element + real(pReal), dimension(:,:), allocatable :: mesh_node, & + mesh_ipVolume ! volume associated with IP + real(pReal), dimension(:,:,:), allocatable :: mesh_ipArea, & ! area of interface to neighboring IP + mesh_ipCenterOfGravity ! center of gravity of IP + real(pReal), dimension(:,:,:,:), allocatable :: mesh_ipAreaNormal ! area normal of interface to neighboring IP + + integer(pInt), dimension(:,:,:), allocatable :: FE_nodesAtIP + integer(pInt), dimension(:,:,:), allocatable :: FE_ipNeighbor + integer(pInt), dimension(:,:,:), allocatable :: FE_subNodeParent + integer(pInt), dimension(:,:,:,:), allocatable :: FE_subNodeOnIPFace + + logical :: noPart ! for cases where the ABAQUS input file does not use part/assembly information + logical, dimension(3) :: mesh_periodicSurface ! flag indicating periodic outer surfaces (used for fluxes) + + integer(pInt) :: hypoelasticTableStyle + integer(pInt) :: initialcondTableStyle + integer(pInt), parameter :: FE_Nelemtypes = 9 + integer(pInt), parameter :: FE_maxNnodes = 8 + integer(pInt), parameter :: FE_maxNsubNodes = 56 + integer(pInt), parameter :: FE_maxNips = 27 + integer(pInt), parameter :: FE_maxNipNeighbors = 6 + integer(pInt), parameter :: FE_maxmaxNnodesAtIP = 8 + integer(pInt), parameter :: FE_NipFaceNodes = 4 + integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_Nnodes = & + (/8, & ! element 7 + 4, & ! element 134 + 4, & ! element 11 + 4, & ! element 27 + 4, & ! element 157 + 6, & ! element 136 + 8, & ! element 21 + 8, & ! element 117 + 8 & ! element 57 (c3d20r == c3d8 --> copy of 7) + /) + integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_NoriginalNodes = & + (/8, & ! element 7 + 4, & ! element 134 + 4, & ! element 11 + 8, & ! element 27 + 4, & ! element 157 + 6, & ! element 136 + 20,& ! element 21 + 8, & ! element 117 + 20 & ! element 57 (c3d20r == c3d8 --> copy of 7) + /) + integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_Nips = & + (/8, & ! element 7 + 1, & ! element 134 + 4, & ! element 11 + 9, & ! element 27 + 4, & ! element 157 + 6, & ! element 136 + 27,& ! element 21 + 1, & ! element 117 + 8 & ! element 57 (c3d20r == c3d8 --> copy of 7) + /) + integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_NipNeighbors = & + (/6, & ! element 7 + 4, & ! element 134 + 4, & ! element 11 + 4, & ! element 27 + 6, & ! element 157 + 6, & ! element 136 + 6, & ! element 21 + 6, & ! element 117 + 6 & ! element 57 (c3d20r == c3d8 --> copy of 7) + /) + integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_NsubNodes = & + (/19,& ! element 7 + 0, & ! element 134 + 5, & ! element 11 + 12,& ! element 27 + 0, & ! element 157 + 15,& ! element 136 + 56,& ! element 21 + 0, & ! element 117 + 19 & ! element 57 (c3d20r == c3d8 --> copy of 7) + /) + integer(pInt), dimension(FE_maxNipNeighbors,FE_Nelemtypes), parameter :: FE_NfaceNodes = & + reshape((/& + 4,4,4,4,4,4, & ! element 7 + 3,3,3,3,0,0, & ! element 134 + 2,2,2,2,0,0, & ! element 11 + 2,2,2,2,0,0, & ! element 27 + 3,3,3,3,0,0, & ! element 157 + 3,4,4,4,3,0, & ! element 136 + 4,4,4,4,4,4, & ! element 21 + 4,4,4,4,4,4, & ! element 117 + 4,4,4,4,4,4 & ! element 57 (c3d20r == c3d8 --> copy of 7) + /),(/FE_maxNipNeighbors,FE_Nelemtypes/)) + integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_maxNnodesAtIP = & + (/1, & ! element 7 + 4, & ! element 134 + 1, & ! element 11 + 2, & ! element 27 + 1, & ! element 157 + 1, & ! element 136 + 4, & ! element 21 + 8, & ! element 117 + 1 & ! element 57 (c3d20r == c3d8 --> copy of 7) + /) + integer(pInt), dimension(FE_NipFaceNodes,FE_maxNipNeighbors,FE_Nelemtypes), parameter :: FE_nodeOnFace = & + reshape((/& + 1,2,3,4 , & ! element 7 + 2,1,5,6 , & + 3,2,6,7 , & + 4,3,7,8 , & + 4,1,5,8 , & + 8,7,6,5 , & + 1,2,3,0 , & ! element 134 + 1,4,2,0 , & + 2,3,4,0 , & + 1,3,4,0 , & + 0,0,0,0 , & + 0,0,0,0 , & + 1,2,0,0 , & ! element 11 + 2,3,0,0 , & + 3,4,0,0 , & + 4,1,0,0 , & + 0,0,0,0 , & + 0,0,0,0 , & + 1,2,0,0 , & ! element 27 + 2,3,0,0 , & + 3,4,0,0 , & + 4,1,0,0 , & + 0,0,0,0 , & + 0,0,0,0 , & + 1,2,3,0 , & ! element 157 + 1,4,2,0 , & + 2,3,4,0 , & + 1,3,4,0 , & + 0,0,0,0 , & + 0,0,0,0 , & + 1,2,3,0 , & ! element 136 + 1,4,5,2 , & + 2,5,6,3 , & + 1,3,6,4 , & + 4,6,5,0 , & + 0,0,0,0 , & + 1,2,3,4 , & ! element 21 + 2,1,5,6 , & + 3,2,6,7 , & + 4,3,7,8 , & + 4,1,5,8 , & + 8,7,6,5 , & + 1,2,3,4 , & ! element 117 + 2,1,5,6 , & + 3,2,6,7 , & + 4,3,7,8 , & + 4,1,5,8 , & + 8,7,6,5 , & + 1,2,3,4 , & ! element 57 (c3d20r == c3d8 --> copy of 7) + 2,1,5,6 , & + 3,2,6,7 , & + 4,3,7,8 , & + 4,1,5,8 , & + 8,7,6,5 & + /),(/FE_NipFaceNodes,FE_maxNipNeighbors,FE_Nelemtypes/)) + + CONTAINS +! --------------------------- +! subroutine mesh_init() +! function mesh_FEtoCPelement(FEid) +! function mesh_build_ipNeighorhood() +! --------------------------- + + +!*********************************************************** +! initialization +!*********************************************************** + subroutine mesh_init (ip,element) + + use mpie_interface + use prec, only: pInt + use IO, only: IO_error,IO_open_InputFile,IO_abaqus_hasNoPart + use FEsolving, only: parallelExecution, FEsolving_execElem, FEsolving_execIP, calcMode, lastMode + + implicit none + + integer(pInt), parameter :: fileUnit = 222 + integer(pInt) e,element,ip + + write(6,*) + write(6,*) '<<<+- mesh init -+>>>' + write(6,*) '$Id: mesh.f90 765 2011-02-16 16:23:08Z MPIE\c.kords $' + write(6,*) + + call mesh_build_FEdata() ! --- get properties of the different types of elements + + if (IO_open_inputFile(fileUnit)) then ! --- parse info from input file... + + select case (FEsolver) + case ('Spectral') + call mesh_spectral_count_nodesAndElements(fileUnit) + call mesh_spectral_count_cpElements() + call mesh_spectral_map_elements() + call mesh_spectral_map_nodes() + call mesh_spectral_count_cpSizes() + call mesh_spectral_build_nodes(fileUnit) + call mesh_spectral_build_elements(fileUnit) + + case ('Marc') + call mesh_marc_get_tableStyles(fileUnit) + call mesh_marc_count_nodesAndElements(fileUnit) + call mesh_marc_count_elementSets(fileUnit) + call mesh_marc_map_elementSets(fileUnit) + call mesh_marc_count_cpElements(fileUnit) + call mesh_marc_map_elements(fileUnit) + call mesh_marc_map_nodes(fileUnit) + call mesh_marc_build_nodes(fileUnit) + call mesh_marc_count_cpSizes(fileunit) + call mesh_marc_build_elements(fileUnit) + call mesh_marc_get_mpieOptions(fileUnit) + case ('Abaqus') + noPart = IO_abaqus_hasNoPart(fileUnit) + call mesh_abaqus_count_nodesAndElements(fileUnit) + call mesh_abaqus_count_elementSets(fileUnit) + call mesh_abaqus_count_materials(fileUnit) + call mesh_abaqus_map_elementSets(fileUnit) + call mesh_abaqus_map_materials(fileUnit) + call mesh_abaqus_count_cpElements(fileUnit) + call mesh_abaqus_map_elements(fileUnit) + call mesh_abaqus_map_nodes(fileUnit) + call mesh_abaqus_build_nodes(fileUnit) + call mesh_abaqus_count_cpSizes(fileunit) + call mesh_abaqus_build_elements(fileUnit) + end select + close (fileUnit) + + call mesh_build_subNodeCoords() + call mesh_build_ipVolumes() + call mesh_build_ipAreas() + call mesh_build_nodeTwins() + call mesh_build_sharedElems() + call mesh_build_ipNeighborhood() + call mesh_tell_statistics() + + parallelExecution = (parallelExecution .and. (mesh_Nelems == mesh_NcpElems)) ! plus potential killer from non-local constitutive + else + call IO_error(101) ! cannot open input file + endif + + FEsolving_execElem = (/1,mesh_NcpElems/) + allocate(FEsolving_execIP(2,mesh_NcpElems)); FEsolving_execIP = 1_pInt + forall (e = 1:mesh_NcpElems) FEsolving_execIP(2,e) = FE_Nips(mesh_element(2,e)) + + allocate(calcMode(mesh_maxNips,mesh_NcpElems)) + write(6,*) '<<<+- mesh init done -+>>>' + calcMode = .false. ! pretend to have collected what first call is asking (F = I) + calcMode(ip,mesh_FEasCP('elem',element)) = .true. ! first ip,el needs to be already pingponged to "calc" + lastMode = .true. ! and its mode is already known... + endsubroutine + + + +!*********************************************************** +! mapping of FE element types to internal representation +!*********************************************************** + function FE_mapElemtype(what) + + use IO, only: IO_lc + + implicit none + + character(len=*), intent(in) :: what + integer(pInt) FE_mapElemtype + + select case (IO_lc(what)) + case ( '7', & + 'c3d8') + FE_mapElemtype = 1 ! Three-dimensional Arbitrarily Distorted Brick + case ('134', & + 'c3d4') + FE_mapElemtype = 2 ! Three-dimensional Four-node Tetrahedron + case ( '11', & + 'cpe4') + FE_mapElemtype = 3 ! Arbitrary Quadrilateral Plane-strain + case ( '27', & + 'cpe8') + FE_mapElemtype = 4 ! Plane Strain, Eight-node Distorted Quadrilateral + case ('157') + FE_mapElemtype = 5 ! Three-dimensional, Low-order, Tetrahedron, Herrmann Formulations + case ('136', & + 'c3d6') + FE_mapElemtype = 6 ! Three-dimensional Arbitrarily Distorted Pentahedral + case ( '21', & + 'c3d20') + FE_mapElemtype = 7 ! Three-dimensional Arbitrarily Distorted quadratic hexahedral + case ( '117', & + '123', & + 'c3d8r') + FE_mapElemtype = 8 ! Three-dimensional Arbitrarily Distorted linear hexahedral with reduced integration + case ( '57', & + 'c3d20r') + FE_mapElemtype = 9 ! Three-dimensional Arbitrarily Distorted quad hexahedral with reduced integration + case default + FE_mapElemtype = 0 ! unknown element --> should raise an error upstream..! + endselect + + endfunction + + + +!*********************************************************** +! FE to CP id mapping by binary search thru lookup array +! +! valid questions are 'elem', 'node' +!*********************************************************** + function mesh_FEasCP(what,id) + + use prec, only: pInt + use IO, only: IO_lc + implicit none + + character(len=*), intent(in) :: what + integer(pInt), intent(in) :: id + integer(pInt), dimension(:,:), pointer :: lookupMap + integer(pInt) mesh_FEasCP, lower,upper,center + + mesh_FEasCP = 0_pInt + select case(IO_lc(what(1:4))) + case('elem') + lookupMap => mesh_mapFEtoCPelem + case('node') + lookupMap => mesh_mapFEtoCPnode + case default + return + endselect + + lower = 1_pInt + upper = size(lookupMap,2) + + ! check at bounds QUESTION is it valid to extend bounds by 1 and just do binary search w/o init check at bounds? + if (lookupMap(1,lower) == id) then + mesh_FEasCP = lookupMap(2,lower) + return + elseif (lookupMap(1,upper) == id) then + mesh_FEasCP = lookupMap(2,upper) + return + endif + + ! binary search in between bounds + do while (upper-lower > 1) + center = (lower+upper)/2 + if (lookupMap(1,center) < id) then + lower = center + elseif (lookupMap(1,center) > id) then + upper = center + else + mesh_FEasCP = lookupMap(2,center) + exit + endif + enddo + return + + endfunction + + +!*********************************************************** +! find face-matching element of same type +!*********************************************************** +subroutine mesh_faceMatch(elem, face ,matchingElem, matchingFace) + +use prec, only: pInt +implicit none + +!*** output variables +integer(pInt), intent(out) :: matchingElem, & ! matching CP element ID + matchingFace ! matching FE face ID + +!*** input variables +integer(pInt), intent(in) :: face, & ! FE face ID + elem ! FE elem ID + +!*** local variables +integer(pInt), dimension(FE_NfaceNodes(face,mesh_element(2,elem))) :: & + myFaceNodes ! global node ids on my face +integer(pInt) myType, & + candidateType, & + candidateElem, & + candidateFace, & + candidateFaceNode, & + minNsharedElems, & + NsharedElems, & + lonelyNode, & + i, & + n, & + dir ! periodicity direction +integer(pInt), dimension(:), allocatable :: element_seen +logical checkTwins + + +minNsharedElems = mesh_maxNsharedElems + 1_pInt ! init to worst case +matchingFace = 0_pInt +matchingElem = 0_pInt ! intialize to "no match found" +myType = mesh_element(2,elem) ! figure elemType + +do n = 1,FE_NfaceNodes(face,myType) ! loop over nodes on face + myFaceNodes(n) = mesh_FEasCP('node',mesh_element(4+FE_nodeOnFace(n,face,myType),elem)) ! CP id of face node + NsharedElems = mesh_sharedElem(1,myFaceNodes(n)) ! figure # shared elements for this node + if (NsharedElems < minNsharedElems) then + minNsharedElems = NsharedElems ! remember min # shared elems + lonelyNode = n ! remember most lonely node + endif +enddo + +allocate(element_seen(minNsharedElems)) +element_seen = 0_pInt + +checkCandidate: do i = 1,minNsharedElems ! iterate over lonelyNode's shared elements + candidateElem = mesh_sharedElem(1+i,myFaceNodes(lonelyNode)) ! present candidate elem + if (all(element_seen /= candidateElem)) then ! element seen for the first time? + element_seen(i) = candidateElem + candidateType = mesh_element(2,candidateElem) ! figure elemType of candidate +checkCandidateFace: do candidateFace = 1,FE_maxNipNeighbors ! check each face of candidate + if (FE_NfaceNodes(candidateFace,candidateType) /= FE_NfaceNodes(face,myType) & ! incompatible face + .or. (candidateElem == elem .and. candidateFace == face)) then ! this is my face + cycle checkCandidateFace + endif + checkTwins = .false. + do n = 1,FE_NfaceNodes(candidateFace,candidateType) ! loop through nodes on face + candidateFaceNode = mesh_FEasCP('node', mesh_element(4+FE_nodeOnFace(n,candidateFace,candidateType),candidateElem)) + if (all(myFaceNodes /= candidateFaceNode)) then ! candidate node does not match any of my face nodes + checkTwins = .true. ! perhaps the twin nodes do match + exit + endif + enddo + if(checkTwins) then +checkCandidateFaceTwins: do dir = 1,3 + do n = 1,FE_NfaceNodes(candidateFace,candidateType) ! loop through nodes on face + candidateFaceNode = mesh_FEasCP('node', mesh_element(4+FE_nodeOnFace(n,candidateFace,candidateType),candidateElem)) + if (all(myFaceNodes /= mesh_nodeTwins(dir,candidateFaceNode))) then ! node twin does not match either + if (dir == 3) then + cycle checkCandidateFace + else + cycle checkCandidateFaceTwins ! try twins in next dimension + endif + endif + enddo + exit checkCandidateFaceTwins + enddo checkCandidateFaceTwins + endif + matchingFace = candidateFace + matchingElem = candidateElem + exit checkCandidate ! found my matching candidate + enddo checkCandidateFace + endif +enddo checkCandidate + +deallocate(element_seen) + +endsubroutine + + +!******************************************************************** +! get properties of different types of finite elements +! +! assign globals: +! FE_nodesAtIP, FE_ipNeighbor, FE_subNodeParent, FE_subNodeOnIPFace +!******************************************************************** + subroutine mesh_build_FEdata () + + use prec, only: pInt + implicit none + + allocate(FE_nodesAtIP(FE_maxmaxNnodesAtIP,FE_maxNips,FE_Nelemtypes)) ; FE_nodesAtIP = 0_pInt + allocate(FE_ipNeighbor(FE_maxNipNeighbors,FE_maxNips,FE_Nelemtypes)) ; FE_ipNeighbor = 0_pInt + allocate(FE_subNodeParent(FE_maxNips,FE_maxNsubNodes,FE_Nelemtypes)) ; FE_subNodeParent = 0_pInt + allocate(FE_subNodeOnIPFace(FE_NipFaceNodes,FE_maxNipNeighbors,FE_maxNips,FE_Nelemtypes)) ; FE_subNodeOnIPFace = 0_pInt + + ! fill FE_nodesAtIP with data + FE_nodesAtIP(:,:FE_Nips(1),1) = & ! element 7 + reshape((/& + 1, & + 2, & + 4, & + 3, & + 5, & + 6, & + 8, & + 7 & + /),(/FE_maxNnodesAtIP(1),FE_Nips(1)/)) + FE_nodesAtIP(:,:FE_Nips(2),2) = & ! element 134 + reshape((/& + 1,2,3,4 & + /),(/FE_maxNnodesAtIP(2),FE_Nips(2)/)) + FE_nodesAtIP(:,:FE_Nips(3),3) = & ! element 11 + reshape((/& + 1, & + 2, & + 4, & + 3 & + /),(/FE_maxNnodesAtIP(3),FE_Nips(3)/)) + FE_nodesAtIP(:,:FE_Nips(4),4) = & ! element 27 + reshape((/& + 1,0, & + 1,2, & + 2,0, & + 1,4, & + 0,0, & + 2,3, & + 4,0, & + 3,4, & + 3,0 & + /),(/FE_maxNnodesAtIP(4),FE_Nips(4)/)) + FE_nodesAtIP(:,:FE_Nips(5),5) = & ! element 157 + reshape((/& + 1, & + 2, & + 3, & + 4 & + /),(/FE_maxNnodesAtIP(5),FE_Nips(5)/)) + FE_nodesAtIP(:,:FE_Nips(6),6) = & ! element 136 + reshape((/& + 1, & + 2, & + 3, & + 4, & + 5, & + 6 & + /),(/FE_maxNnodesAtIP(6),FE_Nips(6)/)) + FE_nodesAtIP(:,:FE_Nips(7),7) = & ! element 21 + reshape((/& + 1,0, 0,0, & + 1,2, 0,0, & + 2,0, 0,0, & + 1,4, 0,0, & + 1,3, 2,4, & + 2,3, 0,0, & + 4,0, 0,0, & + 3,4, 0,0, & + 3,0, 0,0, & + 1,5, 0,0, & + 1,6, 2,5, & + 2,6, 0,0, & + 1,8, 4,5, & + 0,0, 0,0, & + 2,7, 3,6, & + 4,8, 0,0, & + 3,8, 4,7, & + 3,7, 0,0, & + 5,0, 0,0, & + 5,6, 0,0, & + 6,0, 0,0, & + 5,8, 0,0, & + 5,7, 6,8, & + 6,7, 0,0, & + 8,0, 0,0, & + 7,8, 0,0, & + 7,0, 0,0 & + /),(/FE_maxNnodesAtIP(7),FE_Nips(7)/)) +! FE_nodesAtIP(:,:FE_Nips(8),8) = & ! element 117 (c3d8r --> single IP per element, so no need for this mapping) +! reshape((/& +! 1,2,3,4,5,6,7,8 & +! /),(/FE_maxNnodesAtIP(8),FE_Nips(8)/)) + FE_nodesAtIP(:,:FE_Nips(9),9) = & ! element 57 (c3d20r == c3d8 --> copy of 7) + reshape((/& + 1, & + 2, & + 4, & + 3, & + 5, & + 6, & + 8, & + 7 & + /),(/FE_maxNnodesAtIP(9),FE_Nips(9)/)) + + ! fill FE_ipNeighbor with data + FE_ipNeighbor(:FE_NipNeighbors(1),:FE_Nips(1),1) = & ! element 7 + reshape((/& + 2,-5, 3,-2, 5,-1, & + -3, 1, 4,-2, 6,-1, & + 4,-5,-4, 1, 7,-1, & + -3, 3,-4, 2, 8,-1, & + 6,-5, 7,-2,-6, 1, & + -3, 5, 8,-2,-6, 2, & + 8,-5,-4, 5,-6, 3, & + -3, 7,-4, 6,-6, 4 & + /),(/FE_NipNeighbors(1),FE_Nips(1)/)) + FE_ipNeighbor(:FE_NipNeighbors(2),:FE_Nips(2),2) = & ! element 134 + reshape((/& + -1,-2,-3,-4 & + /),(/FE_NipNeighbors(2),FE_Nips(2)/)) + FE_ipNeighbor(:FE_NipNeighbors(3),:FE_Nips(3),3) = & ! element 11 + reshape((/& + 2,-4, 3,-1, & + -2, 1, 4,-1, & + 4,-4,-3, 1, & + -2, 3,-3, 2 & + /),(/FE_NipNeighbors(3),FE_Nips(3)/)) + FE_ipNeighbor(:FE_NipNeighbors(4),:FE_Nips(4),4) = & ! element 27 + reshape((/& + 2,-4, 4,-1, & + 3, 1, 5,-1, & + -2, 2, 6,-1, & + 5,-4, 7, 1, & + 6, 4, 8, 2, & + -2, 5, 9, 3, & + 8,-4,-3, 4, & + 9, 7,-3, 5, & + -2, 8,-3, 6 & + /),(/FE_NipNeighbors(4),FE_Nips(4)/)) + FE_ipNeighbor(:FE_NipNeighbors(5),:FE_Nips(5),5) = & ! element 157 + reshape((/& + 2,-4, 3,-2, 4,-1, & + 3,-2, 1,-3, 4,-1, & + 1,-3, 2,-4, 4,-1, & + 1,-3, 2,-4, 3,-2 & + /),(/FE_NipNeighbors(5),FE_Nips(5)/)) + FE_ipNeighbor(:FE_NipNeighbors(6),:FE_Nips(6),6) = & ! element 136 + reshape((/& + 2,-4, 3,-2, 4,-1, & + -3, 1, 3,-2, 5,-1, & + 2,-4,-3, 1, 6,-1, & + 5,-4, 6,-2,-5, 1, & + -3, 4, 6,-2,-5, 2, & + 5,-4,-3, 4,-5, 3 & + /),(/FE_NipNeighbors(6),FE_Nips(6)/)) + FE_ipNeighbor(:FE_NipNeighbors(7),:FE_Nips(7),7) = & ! element 21 + reshape((/& + 2,-5, 4,-2,10,-1, & + 3, 1, 5,-2,11,-1, & + -3, 2, 6,-2,12,-1, & + 5,-5, 7, 1,13,-1, & + 6, 4, 8, 2,14,-1, & + -3, 5, 9, 3,15,-1, & + 8,-5,-4, 4,16,-1, & + 9, 7,-4, 5,17,-1, & + -3, 8,-4, 6,18,-1, & + 11,-5,13,-2,19, 1, & + 12,10,14,-2,20, 2, & + -3,11,15,-2,21, 3, & + 14,-5,16,10,22, 4, & + 15,13,17,11,23, 5, & + -3,14,18,12,24, 6, & + 17,-5,-4,13,25, 7, & + 18,16,-4,14,26, 8, & + -3,17,-4,15,27, 9, & + 20,-5,22,-2,-6,10, & + 21,19,23,-2,-6,11, & + -3,20,24,-2,-6,12, & + 23,-5,25,19,-6,13, & + 24,22,26,20,-6,14, & + -3,23,27,21,-6,15, & + 26,-5,-4,22,-6,16, & + 27,25,-4,23,-6,17, & + -3,26,-4,24,-6,18 & + /),(/FE_NipNeighbors(7),FE_Nips(7)/)) +FE_ipNeighbor(:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117 + reshape((/& + -3,-5,-4,-2,-6,-1 & + /),(/FE_NipNeighbors(8),FE_Nips(8)/)) + FE_ipNeighbor(:FE_NipNeighbors(9),:FE_Nips(9),9) = & ! element 57 (c3d20r == c3d8 --> copy of 7) + reshape((/& + 2,-5, 3,-2, 5,-1, & + -3, 1, 4,-2, 6,-1, & + 4,-5,-4, 1, 7,-1, & + -3, 3,-4, 2, 8,-1, & + 6,-5, 7,-2,-6, 1, & + -3, 5, 8,-2,-6, 2, & + 8,-5,-4, 5,-6, 3, & + -3, 7,-4, 6,-6, 4 & + /),(/FE_NipNeighbors(9),FE_Nips(9)/)) + + ! fill FE_subNodeParent with data + FE_subNodeParent(:FE_Nips(1),:FE_NsubNodes(1),1) = & ! element 7 + reshape((/& + 1, 2, 0, 0, 0, 0, 0, 0, & + 2, 3, 0, 0, 0, 0, 0, 0, & + 3, 4, 0, 0, 0, 0, 0, 0, & + 4, 1, 0, 0, 0, 0, 0, 0, & + 1, 5, 0, 0, 0, 0, 0, 0, & + 2, 6, 0, 0, 0, 0, 0, 0, & + 3, 7, 0, 0, 0, 0, 0, 0, & + 4, 8, 0, 0, 0, 0, 0, 0, & + 5, 6, 0, 0, 0, 0, 0, 0, & + 6, 7, 0, 0, 0, 0, 0, 0, & + 7, 8, 0, 0, 0, 0, 0, 0, & + 8, 5, 0, 0, 0, 0, 0, 0, & + 1, 2, 3, 4, 0, 0, 0, 0, & + 1, 2, 6, 5, 0, 0, 0, 0, & + 2, 3, 7, 6, 0, 0, 0, 0, & + 3, 4, 8, 7, 0, 0, 0, 0, & + 1, 4, 8, 5, 0, 0, 0, 0, & + 5, 6, 7, 8, 0, 0, 0, 0, & + 1, 2, 3, 4, 5, 6, 7, 8 & + /),(/FE_Nips(1),FE_NsubNodes(1)/)) +!FE_subNodeParent(:FE_Nips(2),:FE_NsubNodes(2),2) ! element 134 has no subnodes + FE_subNodeParent(:FE_Nips(3),:FE_NsubNodes(3),3) = & ! element 11 + reshape((/& + 1, 2, 0, 0, & + 2, 3, 0, 0, & + 3, 4, 0, 0, & + 4, 1, 0, 0, & + 1, 2, 3, 4 & + /),(/FE_Nips(3),FE_NsubNodes(3)/)) + FE_subNodeParent(:FE_Nips(4),:FE_NsubNodes(4),4) = & ! element 27 + reshape((/& + 1, 1, 2, 0, 0, 0, 0, 0, 0, & + 1, 2, 2, 0, 0, 0, 0, 0, 0, & + 2, 2, 3, 0, 0, 0, 0, 0, 0, & + 2, 3, 3, 0, 0, 0, 0, 0, 0, & + 3, 3, 4, 0, 0, 0, 0, 0, 0, & + 3, 4, 4, 0, 0, 0, 0, 0, 0, & + 4, 4, 1, 0, 0, 0, 0, 0, 0, & + 4, 1, 1, 0, 0, 0, 0, 0, 0, & + 1, 1, 1, 1, 2, 2, 4, 4, 3, & + 2, 2, 2, 2, 1, 1, 3, 3, 4, & + 3, 3, 3, 3, 2, 2, 4, 4, 1, & + 4, 4, 4, 4, 1, 1, 3, 3, 2 & + /),(/FE_Nips(4),FE_NsubNodes(4)/)) + !FE_subNodeParent(:FE_Nips(5),:FE_NsubNodes(5),5) = & ! element 157 + ! reshape((/& + ! *still to be defined* + ! /),(/FE_Nips(5),FE_NsubNodes(5)/)) + FE_subNodeParent(:FE_Nips(6),:FE_NsubNodes(6),6) = & ! element 136 + reshape((/& + 1, 2, 0, 0, 0, 0, & + 2, 3, 0, 0, 0, 0, & + 3, 1, 0, 0, 0, 0, & + 1, 4, 0, 0, 0, 0, & + 2, 5, 0, 0, 0, 0, & + 3, 6, 0, 0, 0, 0, & + 4, 5, 0, 0, 0, 0, & + 5, 6, 0, 0, 0, 0, & + 6, 4, 0, 0, 0, 0, & + 1, 2, 3, 0, 0, 0, & + 1, 2, 4, 5, 0, 0, & + 2, 3, 5, 6, 0, 0, & + 1, 3, 4, 6, 0, 0, & + 4, 5, 6, 0, 0, 0, & + 1, 2, 3, 4, 5, 6 & + /),(/FE_Nips(6),FE_NsubNodes(6)/)) + FE_subNodeParent(:FE_Nips(7),:FE_NsubNodes(7),7) = & ! element 21 + reshape((/& + 1, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 1, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 2, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 2, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 3, 3, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 3, 4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 4, 4, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 4, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 1, 1, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 2, 2, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 3, 3, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 4, 4, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 1, 5, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 2, 6, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 3, 7, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 4, 8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 5, 5, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 5, 6, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 6, 6, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 6, 7, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 7, 7, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 7, 8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 8, 8, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 8, 5, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 1, 1, 1, 1, 2, 2, 4, 4, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 2, 2, 2, 2, 1, 1, 3, 3, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 3, 3, 3, 3, 2, 2, 4, 4, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 4, 4, 4, 4, 1, 1, 3, 3, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 1, 1, 1, 1, 2, 2, 5, 5, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 2, 2, 2, 2, 1, 1, 6, 6, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 2, 2, 2, 2, 3, 3, 6, 6, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 3, 3, 3, 3, 2, 2, 7, 7, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 3, 3, 3, 3, 4, 4, 7, 7, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 4, 4, 4, 4, 3, 3, 8, 8, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 4, 4, 4, 4, 1, 1, 8, 8, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 1, 1, 1, 1, 4, 4, 5, 5, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 5, 5, 5, 5, 1, 1, 6, 6, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 6, 6, 6, 6, 2, 2, 5, 5, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 6, 6, 6, 6, 2, 2, 7, 7, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 7, 7, 7, 7, 3, 3, 6, 6, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 7, 7, 7, 7, 3, 3, 8, 8, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 8, 8, 8, 8, 4, 4, 7, 7, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 8, 8, 8, 8, 4, 4, 5, 5, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 5, 5, 5, 5, 1, 1, 8, 8, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 5, 5, 5, 5, 6, 6, 8, 8, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 6, 6, 6, 6, 5, 5, 7, 7, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 7, 7, 7, 7, 6, 6, 8, 8, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 8, 8, 8, 8, 5, 5, 7, 7, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 4, 4, 4, 4, 5, 5, 5, 5, 3, 3, 6, 6, 8, 8, 7, & + 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 3, 3, 3, 3, 6, 6, 6, 6, 4, 4, 5, 5, 7, 7, 8, & + 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 4, 4, 4, 4, 7, 7, 7, 7, 1, 1, 6, 6, 8, 8, 5, & + 4, 4, 4, 4, 4, 4, 4, 4, 1, 1, 1, 1, 3, 3, 3, 3, 8, 8, 8, 8, 2, 2, 5, 5, 7, 7, 6, & + 5, 5, 5, 5, 5, 5, 5, 5, 1, 1, 1, 1, 6, 6, 6, 6, 8, 8, 8, 8, 2, 2, 4, 4, 7, 7, 3, & + 6, 6, 6, 6, 6, 6, 6, 6, 2, 2, 2, 2, 5, 5, 5, 5, 7, 7, 7, 7, 1, 1, 3, 3, 8, 8, 4, & + 7, 7, 7, 7, 7, 7, 7, 7, 3, 3, 3, 3, 6, 6, 6, 6, 8, 8, 8, 8, 2, 2, 4, 4, 5, 5, 1, & + 8, 8, 8, 8, 8, 8, 8, 8, 4, 4, 4, 4, 5, 5, 5, 5, 7, 7, 7, 7, 1, 1, 3, 3, 6, 6, 2 & + /),(/FE_Nips(7),FE_NsubNodes(7)/)) +!FE_subNodeParent(:FE_Nips(8),:FE_NsubNodes(8),8) ! element 117 has no subnodes + FE_subNodeParent(:FE_Nips(9),:FE_NsubNodes(9),9) = & ! element 57 (c3d20r == c3d8 --> copy of 7) + reshape((/& + 1, 2, 0, 0, 0, 0, 0, 0, & + 2, 3, 0, 0, 0, 0, 0, 0, & + 3, 4, 0, 0, 0, 0, 0, 0, & + 4, 1, 0, 0, 0, 0, 0, 0, & + 1, 5, 0, 0, 0, 0, 0, 0, & + 2, 6, 0, 0, 0, 0, 0, 0, & + 3, 7, 0, 0, 0, 0, 0, 0, & + 4, 8, 0, 0, 0, 0, 0, 0, & + 5, 6, 0, 0, 0, 0, 0, 0, & + 6, 7, 0, 0, 0, 0, 0, 0, & + 7, 8, 0, 0, 0, 0, 0, 0, & + 8, 5, 0, 0, 0, 0, 0, 0, & + 1, 2, 3, 4, 0, 0, 0, 0, & + 1, 2, 6, 5, 0, 0, 0, 0, & + 2, 3, 7, 6, 0, 0, 0, 0, & + 3, 4, 8, 7, 0, 0, 0, 0, & + 1, 4, 8, 5, 0, 0, 0, 0, & + 5, 6, 7, 8, 0, 0, 0, 0, & + 1, 2, 3, 4, 5, 6, 7, 8 & + /),(/FE_Nips(9),FE_NsubNodes(9)/)) + + ! fill FE_subNodeOnIPFace with data + FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(1),:FE_Nips(1),1) = & ! element 7 + reshape((/& + 9,21,27,22, & ! 1 + 1,13,25,12, & + 12,25,27,21, & + 1, 9,22,13, & + 13,22,27,25, & + 1,12,21, 9, & + 2,10,23,14, & ! 2 + 9,22,27,21, & + 10,21,27,23, & + 2,14,22, 9, & + 14,23,27,22, & + 2, 9,21,10, & + 11,24,27,21, & ! 3 + 4,12,25,16, & + 4,16,24,11, & + 12,21,27,25, & + 16,25,27,24, & + 4,11,21,12, & + 3,15,23,10, & ! 4 + 11,21,27,24, & + 3,11,24,15, & + 10,23,27,21, & + 15,24,27,23, & + 3,10,21,11, & + 17,22,27,26, & ! 5 + 5,20,25,13, & + 20,26,27,25, & + 5,13,22,17, & + 5,17,26,20, & + 13,25,27,22, & + 6,14,23,18, & ! 6 + 17,26,27,22, & + 18,23,27,26, & + 6,17,22,14, & + 6,18,26,17, & + 14,22,27,23, & + 19,26,27,24, & ! 7 + 8,16,25,20, & + 8,19,24,16, & + 20,25,27,26, & + 8,20,26,19, & + 16,24,27,25, & + 7,18,23,15, & ! 8 + 19,24,27,26, & + 7,15,24,19, & + 18,26,27,23, & + 7,19,26,18, & + 15,23,27,24 & + /),(/FE_NipFaceNodes,FE_NipNeighbors(1),FE_Nips(1)/)) + FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(2),:FE_Nips(2),2) = & ! element 134 + reshape((/& + 1, 1, 3, 2, & ! 1 + 1, 1, 2, 4, & + 2, 2, 3, 4, & + 1, 1, 4, 3 & + /),(/FE_NipFaceNodes,FE_NipNeighbors(2),FE_Nips(2)/)) + FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(3),:FE_Nips(3),3) = & ! element 11 + reshape((/& + 5, 9, 0, 0, & ! 1 + 1, 8, 0, 0, & + 8, 9, 0, 0, & + 1, 5, 0, 0, & + 2, 6, 0, 0, & ! 2 + 5, 9, 0, 0, & + 6, 9, 0, 0, & + 2, 5, 0, 0, & + 3, 6, 0, 0, & ! 3 + 7, 9, 0, 0, & + 3, 7, 0, 0, & + 6, 9, 0, 0, & + 7, 9, 0, 0, & ! 4 + 4, 8, 0, 0, & + 4, 7, 0, 0, & + 8, 9, 0, 0 & + /),(/FE_NipFaceNodes,FE_NipNeighbors(3),FE_Nips(3)/)) + FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(4),:FE_Nips(4),4) = & ! element 27 + reshape((/& + 9,17, 0, 0, & ! 1 + 1,16, 0, 0, & + 16,17, 0, 0, & + 1, 9, 0, 0, & + 10,18, 0, 0, & ! 2 + 9,17, 0, 0, & + 17,18, 0, 0, & + 9,10, 0, 0, & + 2,11, 0, 0, & ! 3 + 10,18, 0, 0, & + 11,18, 0, 0, & + 2,10, 0, 0, & + 17,20, 0, 0, & ! 4 + 15,16, 0, 0, & + 15,20, 0, 0, & + 16,17, 0, 0, & + 18,19, 0, 0, & ! 5 + 17,20, 0, 0, & + 19,20, 0, 0, & + 17,18, 0, 0, & + 11,12, 0, 0, & ! 6 + 18,19, 0, 0, & + 12,19, 0, 0, & + 11,18, 0, 0, & + 14,20, 0, 0, & ! 7 + 4,15, 0, 0, & + 4,14, 0, 0, & + 15,20, 0, 0, & + 13,19, 0, 0, & ! 8 + 14,20, 0, 0, & + 13,14, 0, 0, & + 19,20, 0, 0, & + 3,12, 0, 0, & ! 9 + 13,19, 0, 0, & + 3,13, 0, 0, & + 12,19, 0, 0 & + /),(/FE_NipFaceNodes,FE_NipNeighbors(4),FE_Nips(4)/)) + !FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(5),:FE_Nips(5),5) = & ! element 157 + ! reshape((/& + ! *still to be defined* + ! /),(/FE_NipFaceNodes,FE_NipNeighbors(5),FE_Nips(5)/)) + FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(6),:FE_Nips(6),6) = & ! element 136 + reshape((/& + 7,16,21,17, & ! 1 + 1,10,19, 9, & + 9,19,21,16, & + 1, 7,17,10, & + 10,17,21,19, & + 1, 9,16, 7, & + 2, 8,18,11, & ! 2 + 7,17,21,16, & + 8,16,21,18, & + 2,11,17, 7, & + 11,18,21,17, & + 2, 7,16, 8, & + 8,18,21,16, & ! 3 + 3, 9,19,12, & + 3,12,18, 8, & + 9,16,21,19, & + 12,19,21,18, & + 3, 8,16, 9, & + 13,17,21,20, & ! 4 + 4,15,19,10, & + 15,20,21,19, & + 4,10,17,13, & + 4,13,20,15, & + 10,19,21,17, & + 5,11,18,14, & ! 5 + 13,20,21,17, & + 14,18,21,20, & + 5,13,17,11, & + 5,14,20,13, & + 11,17,21,18, & + 14,20,21,18, & ! 6 + 6,12,19,15, & + 6,14,18,12, & + 15,19,21,20, & + 6,15,20,14, & + 12,18,21,19 & + /),(/FE_NipFaceNodes,FE_NipNeighbors(6),FE_Nips(6)/)) + FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(7),:FE_Nips(7),7) = & ! element 21 + reshape((/& + 9,33,57,37, & ! 1 + 1,17,44,16, & + 33,16,44,57, & + 1, 9,37,17, & + 17,37,57,44, & + 1,16,33, 9, & + 10,34,58,38, & ! 2 + 9,37,57,33, & + 34,33,57,58, & + 9,10,38,37, & + 37,38,58,57, & + 9,33,34,10, & + 2,11,39,18, & ! 3 + 10,38,58,34, & + 11,34,58,39, & + 10, 2,18,38, & + 38,18,39,58, & + 10,34,11, 2, & + 33,36,60,57, & ! 4 + 16,44,43,15, & + 36,15,43,60, & + 16,33,57,44, & + 44,57,60,43, & + 16,15,36,33, & + 34,35,59,58, & ! 5 + 33,57,60,36, & + 35,36,60,59, & + 33,34,58,57, & + 57,58,59,60, & + 33,36,35,34, & + 11,12,40,39, & ! 6 + 34,58,59,35, & + 12,35,59,40, & + 34,11,39,58, & + 58,39,40,59, & + 34,35,12,11, & + 36,14,42,60, & ! 7 + 15,43,20, 4, & + 14, 4,20,42, & + 15,36,60,43, & + 43,60,42,20, & + 15, 4,14,36, & + 35,13,41,59, & ! 8 + 36,60,42,14, & + 13,14,42,41, & + 36,35,59,60, & + 60,59,41,42, & + 36,14,13,35, & + 12, 3,19,40, & ! 9 + 35,59,41,13, & + 3,13,41,19, & + 35,12,40,59, & + 59,40,19,41, & + 35,13, 3,12, & + 37,57,61,45, & ! 10 + 17,21,52,44, & + 57,44,52,61, & + 17,37,45,21, & + 21,45,61,52, & + 17,44,57,37, & + 38,58,62,46, & ! 11 + 37,45,61,57, & + 58,57,61,62, & + 37,38,46,45, & + 45,46,62,61, & + 37,57,58,38, & + 18,39,47,22, & ! 12 + 38,46,62,58, & + 39,58,62,47, & + 38,18,22,46, & + 46,22,47,62, & + 38,58,39,18, & + 57,60,64,61, & ! 13 + 44,52,51,43, & + 60,43,51,64, & + 44,57,61,52, & + 52,61,64,51, & + 44,43,60,57, & + 58,59,63,62, & ! 14 + 57,61,64,60, & + 59,60,64,63, & + 57,58,62,61, & + 61,62,63,64, & + 57,60,59,58, & + 39,40,48,47, & ! 15 + 58,62,63,59, & + 40,59,63,48, & + 58,39,47,62, & + 62,47,48,63, & + 58,59,40,39, & + 60,42,50,64, & ! 16 + 43,51,24,20, & + 42,20,24,50, & + 43,60,64,51, & + 51,64,50,24, & + 43,20,42,60, & + 59,41,49,63, & ! 17 + 60,64,50,42, & + 41,42,50,49, & + 60,59,63,64, & + 64,63,49,50, & + 60,42,41,59, & + 40,19,23,48, & ! 18 + 59,63,49,41, & + 19,41,49,23, & + 59,40,48,63, & + 63,48,23,49, & + 59,41,19,40, & + 45,61,53,25, & ! 19 + 21, 5,32,52, & + 61,52,32,53, & + 21,45,25, 5, & + 5,25,53,32, & + 21,52,61,45, & + 46,62,54,26, & ! 20 + 45,25,53,61, & + 62,61,53,54, & + 45,46,26,25, & + 25,26,54,53, & + 45,61,62,46, & + 22,47,27, 6, & ! 21 + 46,26,54,62, & + 47,62,54,27, & + 46,22, 6,26, & + 26, 6,27,54, & + 46,62,47,22, & + 61,64,56,53, & ! 22 + 52,32,31,51, & + 64,51,31,56, & + 52,61,53,32, & + 32,53,56,31, & + 52,51,64,61, & + 62,63,55,54, & ! 23 + 61,53,56,64, & + 63,64,56,55, & + 61,62,54,53, & + 53,54,55,56, & + 61,64,63,62, & + 47,48,28,27, & ! 24 + 62,54,55,63, & + 48,63,55,28, & + 62,47,27,54, & + 54,27,28,55, & + 62,63,48,47, & + 64,50,30,56, & ! 25 + 51,31, 8,24, & + 50,24, 8,30, & + 51,64,56,31, & + 31,56,30, 8, & + 51,24,50,64, & + 63,49,29,55, & ! 26 + 64,56,30,50, & + 49,50,30,29, & + 64,63,55,56, & + 56,55,29,30, & + 64,50,49,63, & + 48,23, 7,28, & ! 27 + 63,55,29,49, & + 23,49,29, 7, & + 63,48,28,55, & + 55,28, 7,29, & + 63,49,23,48 & + /),(/FE_NipFaceNodes,FE_NipNeighbors(7),FE_Nips(7)/)) + FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117 + reshape((/& + 2, 3, 7, 6, & ! 1 + 1, 5, 8, 4, & + 3, 4, 8, 7, & + 1, 2, 6, 5, & + 5, 6, 7, 8, & + 1, 4, 3, 2 & + /),(/FE_NipFaceNodes,FE_NipNeighbors(8),FE_Nips(8)/)) + FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(9),:FE_Nips(9),9) = & ! element 57 (c3d20r == c3d8 --> copy of 7) + reshape((/& + 9,21,27,22, & ! 1 + 1,13,25,12, & + 12,25,27,21, & + 1, 9,22,13, & + 13,22,27,25, & + 1,12,21, 9, & + 2,10,23,14, & ! 2 + 9,22,27,21, & + 10,21,27,23, & + 2,14,22, 9, & + 14,23,27,22, & + 2, 9,21,10, & + 11,24,27,21, & ! 3 + 4,12,25,16, & + 4,16,24,11, & + 12,21,27,25, & + 16,25,27,24, & + 4,11,21,12, & + 3,15,23,10, & ! 4 + 11,21,27,24, & + 3,11,24,15, & + 10,23,27,21, & + 15,24,27,23, & + 3,10,21,11, & + 17,22,27,26, & ! 5 + 5,20,25,13, & + 20,26,27,25, & + 5,13,22,17, & + 5,17,26,20, & + 13,25,27,22, & + 6,14,23,18, & ! 6 + 17,26,27,22, & + 18,23,27,26, & + 6,17,22,14, & + 6,18,26,17, & + 14,22,27,23, & + 19,26,27,24, & ! 7 + 8,16,25,20, & + 8,19,24,16, & + 20,25,27,26, & + 8,20,26,19, & + 16,24,27,25, & + 7,18,23,15, & ! 8 + 19,24,27,26, & + 7,15,24,19, & + 18,26,27,23, & + 7,19,26,18, & + 15,23,27,24 & + /),(/FE_NipFaceNodes,FE_NipNeighbors(9),FE_Nips(9)/)) + + return + + endsubroutine + + +!******************************************************************** +! figure out table styles (Marc only) +! +! initialcondTableStyle, hypoelasticTableStyle +!******************************************************************** + subroutine mesh_marc_get_tableStyles (unit) + + use prec, only: pInt + use IO + implicit none + + integer(pInt), parameter :: maxNchunks = 6 + integer(pInt), dimension (1+2*maxNchunks) :: pos + + integer(pInt) unit + character(len=300) line + + initialcondTableStyle = 0_pInt + hypoelasticTableStyle = 0_pInt + +610 FORMAT(A300) + + rewind(unit) + do + read (unit,610,END=620) line + pos = IO_stringPos(line,maxNchunks) + + if ( IO_lc(IO_stringValue(line,pos,1)) == 'table' .and. pos(1) .GT. 5) then + initialcondTableStyle = IO_intValue(line,pos,4) + hypoelasticTableStyle = IO_intValue(line,pos,5) + exit + endif + enddo + +620 return + + endsubroutine + + +!******************************************************************** +! get any additional mpie options from input file (Marc only) +! +! mesh_periodicSurface +!******************************************************************** +subroutine mesh_marc_get_mpieOptions(unit) + +use prec, only: pInt +use IO +implicit none + +integer(pInt), intent(in) :: unit + +integer(pInt), parameter :: maxNchunks = 5 +integer(pInt), dimension (1+2*maxNchunks) :: pos +integer(pInt) chunk +character(len=300) line + +mesh_periodicSurface = (/.false., .false., .false./) + +610 FORMAT(A300) + +rewind(unit) +do + read (unit,610,END=620) line + pos = IO_stringPos(line,maxNchunks) + + if (IO_lc(IO_stringValue(line,pos,1)) == '$mpie') then ! found keyword for user defined input + if (IO_lc(IO_stringValue(line,pos,2)) == 'periodic' & ! found keyword 'periodic' + .and. pos(1) > 3) then ! and there is at least one more chunk to read + do chunk = 2,pos(1) ! loop through chunks (skipping the keyword) + select case(IO_stringValue(line,pos,chunk)) ! chunk matches keyvalues x,y or z? + case('x') + mesh_periodicSurface(1) = .true. + case('y') + mesh_periodicSurface(2) = .true. + case('z') + mesh_periodicSurface(3) = .true. + end select + enddo + endif + endif +enddo + +620 return + +endsubroutine + + +!******************************************************************** +! count overall number of nodes and elements in mesh +! +! mesh_Nelems, mesh_Nnodes +!******************************************************************** + subroutine mesh_spectral_count_nodesAndElements (unit) + + use prec, only: pInt + use IO + implicit none + + integer(pInt), parameter :: maxNchunks = 7 + integer(pInt), dimension (1+2*maxNchunks) :: pos + integer(pInt) a,b,c,i + + integer(pInt) unit + character(len=1024) line + + mesh_Nnodes = 0_pInt + mesh_Nelems = 0_pInt + + rewind(unit) + do + read(unit,'(a1024)',END=100) line + if (IO_isBlank(line)) cycle ! skip empty lines + pos = IO_stringPos(line,maxNchunks) + + if ( IO_lc(IO_StringValue(line,pos,1)) == 'resolution') then + do i = 2,6,2 + select case (IO_lc(IO_stringValue(line,pos,i))) + case('a') + a = IO_intValue(line,pos,i+1) + case('b') + b = IO_intValue(line,pos,i+1) + case('c') + c = IO_intValue(line,pos,i+1) + end select + enddo + mesh_Nelems = a * b * c + mesh_Nnodes = (1 + a)*(1 + b)*(1 + c) + exit + endif + enddo + +100 return + + endsubroutine + +!******************************************************************** +! count overall number of nodes and elements in mesh +! +! mesh_Nelems, mesh_Nnodes +!******************************************************************** + subroutine mesh_marc_count_nodesAndElements (unit) + + use prec, only: pInt + use IO + implicit none + + integer(pInt), parameter :: maxNchunks = 4 + integer(pInt), dimension (1+2*maxNchunks) :: pos + + integer(pInt) unit + character(len=300) line + + mesh_Nnodes = 0_pInt + mesh_Nelems = 0_pInt + +610 FORMAT(A300) + + rewind(unit) + do + read (unit,610,END=620) line + pos = IO_stringPos(line,maxNchunks) + + if ( IO_lc(IO_StringValue(line,pos,1)) == 'sizing') then + mesh_Nelems = IO_IntValue (line,pos,3) + mesh_Nnodes = IO_IntValue (line,pos,4) + exit + endif + enddo + +620 return + + endsubroutine + +!******************************************************************** +! count overall number of nodes and elements in mesh +! +! mesh_Nelems, mesh_Nnodes +!******************************************************************** + subroutine mesh_abaqus_count_nodesAndElements (unit) + + use prec, only: pInt + use IO + implicit none + + integer(pInt), parameter :: maxNchunks = 2 + integer(pInt), dimension (1+2*maxNchunks) :: pos + character(len=300) line + + integer(pInt) unit + logical inPart + + mesh_Nnodes = 0 + mesh_Nelems = 0 + +610 FORMAT(A300) + + inPart = .false. + rewind(unit) + do + read (unit,610,END=620) line + pos = IO_stringPos(line,maxNchunks) + if ( IO_lc(IO_stringValue(line,pos,1)) == '*part' ) inPart = .true. + if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. & + IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false. + + if (inPart .or. noPart) then + select case ( IO_lc(IO_stringValue(line,pos,1))) + case('*node') + if( & + IO_lc(IO_stringValue(line,pos,2)) /= 'output' .and. & + IO_lc(IO_stringValue(line,pos,2)) /= 'print' .and. & + IO_lc(IO_stringValue(line,pos,2)) /= 'file' .and. & + IO_lc(IO_stringValue(line,pos,2)) /= 'response' & + ) & + mesh_Nnodes = mesh_Nnodes + IO_countDataLines(unit) + case('*element') + if( & + IO_lc(IO_stringValue(line,pos,2)) /= 'output' .and. & + IO_lc(IO_stringValue(line,pos,2)) /= 'matrix' .and. & + IO_lc(IO_stringValue(line,pos,2)) /= 'response' & + ) then + mesh_Nelems = mesh_Nelems + IO_countDataLines(unit) + endif + endselect + endif + enddo + +620 if (mesh_Nnodes < 2) call IO_error(900) + if (mesh_Nelems == 0) call IO_error(901) + + return + + endsubroutine + + +!******************************************************************** +! count overall number of element sets in mesh +! +! mesh_NelemSets, mesh_maxNelemInSet +!******************************************************************** + subroutine mesh_marc_count_elementSets (unit) + + use prec, only: pInt + use IO + implicit none + + integer(pInt), parameter :: maxNchunks = 2 + integer(pInt), dimension (1+2*maxNchunks) :: pos + + integer(pInt) unit + character(len=300) line + + mesh_NelemSets = 0_pInt + mesh_maxNelemInSet = 0_pInt + +610 FORMAT(A300) + + rewind(unit) + do + read (unit,610,END=620) line + pos = IO_stringPos(line,maxNchunks) + + if ( IO_lc(IO_StringValue(line,pos,1)) == 'define' .and. & + IO_lc(IO_StringValue(line,pos,2)) == 'element' ) then + mesh_NelemSets = mesh_NelemSets + 1_pInt + mesh_maxNelemInSet = max(mesh_maxNelemInSet, & + IO_countContinousIntValues(unit)) + endif + enddo + +620 return + + endsubroutine + + +!******************************************************************** +! count overall number of element sets in mesh +! +! mesh_NelemSets, mesh_maxNelemInSet +!******************************************************************** + subroutine mesh_abaqus_count_elementSets (unit) + + use prec, only: pInt + use IO + implicit none + + integer(pInt), parameter :: maxNchunks = 2 + integer(pInt), dimension (1+2*maxNchunks) :: pos + character(len=300) line + + integer(pInt) unit + logical inPart + + mesh_NelemSets = 0_pInt + mesh_maxNelemInSet = mesh_Nelems ! have to be conservative, since Abaqus allows for recursive definitons + +610 FORMAT(A300) + + inPart = .false. + rewind(unit) + do + read (unit,610,END=620) line + pos = IO_stringPos(line,maxNchunks) + if ( IO_lc(IO_stringValue(line,pos,1)) == '*part' ) inPart = .true. + if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. & + IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false. + + if ( (inPart .or. noPart) .and. IO_lc(IO_stringValue(line,pos,1)) == '*elset' ) & + mesh_NelemSets = mesh_NelemSets + 1_pInt + enddo + +620 continue + if (mesh_NelemSets == 0) call IO_error(902) + + return + endsubroutine + + +!******************************************************************** +! count overall number of solid sections sets in mesh (Abaqus only) +! +! mesh_Nmaterials +!******************************************************************** + subroutine mesh_abaqus_count_materials (unit) + + use prec, only: pInt + use IO + implicit none + + integer(pInt), parameter :: maxNchunks = 2 + integer(pInt), dimension (1+2*maxNchunks) :: pos + character(len=300) line + + integer(pInt) unit + logical inPart + + mesh_Nmaterials = 0_pInt + +610 FORMAT(A300) + + inPart = .false. + rewind(unit) + do + read (unit,610,END=620) line + pos = IO_stringPos(line,maxNchunks) + if ( IO_lc(IO_stringValue(line,pos,1)) == '*part' ) inPart = .true. + if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. & + IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false. + + if ( (inPart .or. noPart) .and. & + IO_lc(IO_StringValue(line,pos,1)) == '*solid' .and. & + IO_lc(IO_StringValue(line,pos,2)) == 'section' ) & + mesh_Nmaterials = mesh_Nmaterials + 1_pInt + enddo + +620 if (mesh_Nmaterials == 0) call IO_error(903) + + return + endsubroutine + + +!******************************************************************** +! count overall number of cpElements in mesh +! +! mesh_NcpElems +!******************************************************************** + subroutine mesh_spectral_count_cpElements () + + use prec, only: pInt + implicit none + + mesh_NcpElems = mesh_Nelems + return + + endsubroutine + + +!******************************************************************** +! count overall number of cpElements in mesh +! +! mesh_NcpElems +!******************************************************************** + subroutine mesh_marc_count_cpElements (unit) + + use prec, only: pInt + use IO + implicit none + + integer(pInt), parameter :: maxNchunks = 1 + integer(pInt), dimension (1+2*maxNchunks) :: pos + + integer(pInt) unit,i + character(len=300) line + + mesh_NcpElems = 0_pInt + +610 FORMAT(A300) + + rewind(unit) + do + read (unit,610,END=620) line + pos = IO_stringPos(line,maxNchunks) + + if ( IO_lc(IO_stringValue(line,pos,1)) == 'hypoelastic') then + do i=1,3+hypoelasticTableStyle ! Skip 3 or 4 lines + read (unit,610,END=620) line + enddo + mesh_NcpElems = mesh_NcpElems + IO_countContinousIntValues(unit) + exit + endif + enddo + +620 return + + endsubroutine + + +!******************************************************************** +! count overall number of cpElements in mesh +! +! mesh_NcpElems +!******************************************************************** + subroutine mesh_abaqus_count_cpElements (unit) + + use prec, only: pInt + use IO + implicit none + + integer(pInt), parameter :: maxNchunks = 2 + integer(pInt), dimension (1+2*maxNchunks) :: pos + character(len=300) line + + integer(pInt) unit,i + logical materialFound + character (len=64) materialName,elemSetName + + mesh_NcpElems = 0 + +610 FORMAT(A300) + + rewind(unit) + do + read (unit,610,END=620) line + pos = IO_stringPos(line,maxNchunks) + select case ( IO_lc(IO_stringValue(line,pos,1)) ) + case('*material') + materialName = IO_extractValue(IO_lc(IO_stringValue(line,pos,2)),'name') ! extract name=value + materialFound = materialName /= '' ! valid name? + case('*user') + if (IO_lc(IO_StringValue(line,pos,2)) == 'material' .and. materialFound) then + do i = 1,mesh_Nmaterials ! look thru material names + if (materialName == mesh_nameMaterial(i)) exit ! found one + enddo + if (i <= mesh_Nmaterials) then ! found one? + elemSetName = mesh_mapMaterial(i) ! take corresponding elemSet + do i = 1,mesh_NelemSets ! look thru all elemSet definitions + if (elemSetName == mesh_nameElemSet(i)) & ! matched? + mesh_NcpElems = mesh_NcpElems + mesh_mapElemSet(1,i) ! add those elem count + enddo + endif + materialFound = .false. + endif + endselect + enddo + +620 if (mesh_NcpElems == 0) call IO_error(906) + + return + endsubroutine + + +!******************************************************************** +! map element sets +! +! allocate globals: mesh_nameElemSet, mesh_mapElemSet +!******************************************************************** + subroutine mesh_marc_map_elementSets (unit) + + use prec, only: pInt + use IO + + implicit none + + integer(pInt), parameter :: maxNchunks = 4 + integer(pInt), dimension (1+2*maxNchunks) :: pos + character(len=300) line + + integer(pInt) unit,elemSet + + allocate (mesh_nameElemSet(mesh_NelemSets)) ; mesh_nameElemSet = '' + allocate (mesh_mapElemSet(1+mesh_maxNelemInSet,mesh_NelemSets)) ; mesh_mapElemSet = 0_pInt + +610 FORMAT(A300) + + elemSet = 0_pInt + rewind(unit) + do + read (unit,610,END=640) line + pos = IO_stringPos(line,maxNchunks) + if( (IO_lc(IO_stringValue(line,pos,1)) == 'define' ) .and. & + (IO_lc(IO_stringValue(line,pos,2)) == 'element' ) ) then + elemSet = elemSet+1 + mesh_nameElemSet(elemSet) = IO_stringValue(line,pos,4) + mesh_mapElemSet(:,elemSet) = IO_continousIntValues(unit,mesh_maxNelemInSet,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets) + endif + enddo + +640 return + + endsubroutine + + +!******************************************************************** +! Build element set mapping +! +! allocate globals: mesh_nameElemSet, mesh_mapElemSet +!******************************************************************** + subroutine mesh_abaqus_map_elementSets (unit) + + use prec, only: pInt + use IO + + implicit none + + integer(pInt), parameter :: maxNchunks = 4 + integer(pInt), dimension (1+2*maxNchunks) :: pos + character(len=300) line + + integer(pInt) unit,elemSet,i + logical inPart + + allocate (mesh_nameElemSet(mesh_NelemSets)) ; mesh_nameElemSet = '' + allocate (mesh_mapElemSet(1+mesh_maxNelemInSet,mesh_NelemSets)) ; mesh_mapElemSet = 0_pInt + +610 FORMAT(A300) + + elemSet = 0_pInt + inPart = .false. + rewind(unit) + do + read (unit,610,END=640) line + pos = IO_stringPos(line,maxNchunks) + if ( IO_lc(IO_stringValue(line,pos,1)) == '*part' ) inPart = .true. + if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. & + IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false. + + if ( (inPart .or. noPart) .and. IO_lc(IO_stringValue(line,pos,1)) == '*elset' ) then + elemSet = elemSet + 1_pInt + mesh_nameElemSet(elemSet) = IO_extractValue(IO_lc(IO_stringValue(line,pos,2)),'elset') + mesh_mapElemSet(:,elemSet) = IO_continousIntValues(unit,mesh_Nelems,mesh_nameElemSet,mesh_mapElemSet,elemSet-1) + endif + enddo + +640 do i = 1,elemSet +! write(6,*)'elemSetName: ',mesh_nameElemSet(i) +! write(6,*)'elems in Elset',mesh_mapElemSet(:,i) + if (mesh_mapElemSet(1,i) == 0) call IO_error(ID=904,ext_msg=mesh_nameElemSet(i)) + enddo + + return + endsubroutine + + +!******************************************************************** +! map solid section (Abaqus only) +! +! allocate globals: mesh_nameMaterial, mesh_mapMaterial +!******************************************************************** + subroutine mesh_abaqus_map_materials (unit) + + use prec, only: pInt + use IO + implicit none + + integer(pInt), parameter :: maxNchunks = 20 + integer(pInt), dimension (1+2*maxNchunks) :: pos + character(len=300) line + + integer(pInt) unit,i,count + logical inPart,materialFound + character(len=64) elemSetName,materialName + + allocate (mesh_nameMaterial(mesh_Nmaterials)) ; mesh_nameMaterial = '' + allocate (mesh_mapMaterial(mesh_Nmaterials)) ; mesh_mapMaterial = '' + +610 FORMAT(A300) + + count = 0 + inPart = .false. + rewind(unit) + do + read (unit,610,END=620) line + pos = IO_stringPos(line,maxNchunks) + if ( IO_lc(IO_stringValue(line,pos,1)) == '*part' ) inPart = .true. + if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. & + IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false. + + if ( (inPart .or. noPart) .and. & + IO_lc(IO_StringValue(line,pos,1)) == '*solid' .and. & + IO_lc(IO_StringValue(line,pos,2)) == 'section' ) then + + elemSetName = '' + materialName = '' + + do i = 3,pos(1) + if (IO_extractValue(IO_lc(IO_stringValue(line,pos,i)),'elset') /= '') & + elemSetName = IO_extractValue(IO_lc(IO_stringValue(line,pos,i)),'elset') + if (IO_extractValue(IO_lc(IO_stringValue(line,pos,i)),'material') /= '') & + materialName = IO_extractValue(IO_lc(IO_stringValue(line,pos,i)),'material') + enddo + + if (elemSetName /= '' .and. materialName /= '') then + count = count + 1_pInt + mesh_nameMaterial(count) = materialName ! name of material used for this section + mesh_mapMaterial(count) = elemSetName ! mapped to respective element set + endif + endif + enddo + +620 if (count==0) call IO_error(905) + do i=1,count +! write(6,*)'name of materials: ',i,mesh_nameMaterial(i) +! write(6,*)'name of elemSets: ',i,mesh_mapMaterial(i) + if (mesh_nameMaterial(i)=='' .or. mesh_mapMaterial(i)=='') call IO_error(905) + enddo + + return + endsubroutine + + + +!******************************************************************** +! map nodes from FE id to internal (consecutive) representation +! +! allocate globals: mesh_mapFEtoCPnode +!******************************************************************** + subroutine mesh_spectral_map_nodes () + + use prec, only: pInt + + implicit none + integer(pInt) i + + allocate (mesh_mapFEtoCPnode(2,mesh_Nnodes)) ; mesh_mapFEtoCPnode = 0_pInt + + forall (i = 1:mesh_Nnodes) & + mesh_mapFEtoCPnode(:,i) = i + + return + + endsubroutine + + + +!******************************************************************** +! map nodes from FE id to internal (consecutive) representation +! +! allocate globals: mesh_mapFEtoCPnode +!******************************************************************** + subroutine mesh_marc_map_nodes (unit) + + use prec, only: pInt + use math, only: qsort + use IO + + implicit none + + integer(pInt), parameter :: maxNchunks = 1 + integer(pInt), dimension (1+2*maxNchunks) :: pos + character(len=300) line + + integer(pInt), dimension (mesh_Nnodes) :: node_count + integer(pInt) unit,i + + allocate (mesh_mapFEtoCPnode(2,mesh_Nnodes)) ; mesh_mapFEtoCPnode = 0_pInt + +610 FORMAT(A300) + + node_count = 0_pInt + + rewind(unit) + do + read (unit,610,END=650) line + pos = IO_stringPos(line,maxNchunks) + if( IO_lc(IO_stringValue(line,pos,1)) == 'coordinates' ) then + read (unit,610,END=650) line ! skip crap line + do i = 1,mesh_Nnodes + read (unit,610,END=650) line + mesh_mapFEtoCPnode(1,i) = IO_fixedIntValue (line,(/0,10/),1) + mesh_mapFEtoCPnode(2,i) = i + enddo + exit + endif + enddo + +650 call qsort(mesh_mapFEtoCPnode,1,size(mesh_mapFEtoCPnode,2)) + + return + + endsubroutine + + + +!******************************************************************** +! map nodes from FE id to internal (consecutive) representation +! +! allocate globals: mesh_mapFEtoCPnode +!******************************************************************** + subroutine mesh_abaqus_map_nodes (unit) + + use prec, only: pInt + use math, only: qsort + use IO + + implicit none + + integer(pInt), parameter :: maxNchunks = 2 + integer(pInt), dimension (1+2*maxNchunks) :: pos + character(len=300) line + + integer(pInt) unit,i,count,cpNode + logical inPart + + allocate (mesh_mapFEtoCPnode(2,mesh_Nnodes)) ; mesh_mapFEtoCPnode = 0_pInt + +610 FORMAT(A300) + + cpNode = 0_pInt + inPart = .false. + rewind(unit) + do + read (unit,610,END=650) line + pos = IO_stringPos(line,maxNchunks) + if ( IO_lc(IO_stringValue(line,pos,1)) == '*part' ) inPart = .true. + if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. & + IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false. + + if( (inPart .or. noPart) .and. & + IO_lc(IO_stringValue(line,pos,1)) == '*node' .and. & + ( IO_lc(IO_stringValue(line,pos,2)) /= 'output' .and. & + IO_lc(IO_stringValue(line,pos,2)) /= 'print' .and. & + IO_lc(IO_stringValue(line,pos,2)) /= 'file' .and. & + IO_lc(IO_stringValue(line,pos,2)) /= 'response' ) & + ) then + count = IO_countDataLines(unit) + do i = 1,count + backspace(unit) + enddo + do i = 1,count + read (unit,610,END=650) line + pos = IO_stringPos(line,maxNchunks) + cpNode = cpNode + 1_pInt + mesh_mapFEtoCPnode(1,cpNode) = IO_intValue(line,pos,1) + mesh_mapFEtoCPnode(2,cpNode) = cpNode + enddo + endif + enddo + +650 call qsort(mesh_mapFEtoCPnode,1,size(mesh_mapFEtoCPnode,2)) + + if (size(mesh_mapFEtoCPnode) == 0) call IO_error(908) + return + + endsubroutine + + +!******************************************************************** +! map elements from FE id to internal (consecutive) representation +! +! allocate globals: mesh_mapFEtoCPelem +!******************************************************************** + subroutine mesh_spectral_map_elements () + + use prec, only: pInt + + implicit none + integer(pInt) i + + allocate (mesh_mapFEtoCPelem(2,mesh_NcpElems)) ; mesh_mapFEtoCPelem = 0_pInt + + forall (i = 1:mesh_NcpElems) & + mesh_mapFEtoCPelem(:,i) = i + + return + + endsubroutine + + + +!******************************************************************** +! map elements from FE id to internal (consecutive) representation +! +! allocate globals: mesh_mapFEtoCPelem +!******************************************************************** + subroutine mesh_marc_map_elements (unit) + + use prec, only: pInt + use math, only: qsort + use IO + + implicit none + + integer(pInt), parameter :: maxNchunks = 1 + integer(pInt), dimension (1+2*maxNchunks) :: pos + character(len=300) line + + integer(pInt), dimension (1+mesh_NcpElems) :: contInts + integer(pInt) unit,i,cpElem + + allocate (mesh_mapFEtoCPelem(2,mesh_NcpElems)) ; mesh_mapFEtoCPelem = 0_pInt + +610 FORMAT(A300) + + cpElem = 0_pInt + rewind(unit) + do + read (unit,610,END=660) line + pos = IO_stringPos(line,maxNchunks) + if( IO_lc(IO_stringValue(line,pos,1)) == 'hypoelastic' ) then + do i=1,3+hypoelasticTableStyle ! skip three (or four if new table style!) lines + read (unit,610,END=660) line + enddo + contInts = IO_continousIntValues(unit,mesh_NcpElems,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets) + do i = 1,contInts(1) + cpElem = cpElem+1 + mesh_mapFEtoCPelem(1,cpElem) = contInts(1+i) + mesh_mapFEtoCPelem(2,cpElem) = cpElem + enddo + endif + enddo + +660 call qsort(mesh_mapFEtoCPelem,1,size(mesh_mapFEtoCPelem,2)) ! should be mesh_NcpElems + + return + endsubroutine + + +!******************************************************************** +! map elements from FE id to internal (consecutive) representation +! +! allocate globals: mesh_mapFEtoCPelem +!******************************************************************** + subroutine mesh_abaqus_map_elements (unit) + + use prec, only: pInt + use math, only: qsort + use IO + + implicit none + + integer(pInt), parameter :: maxNchunks = 2 + integer(pInt), dimension (1+2*maxNchunks) :: pos + character(len=300) line + + integer(pInt) unit,i,j,cpElem + logical materialFound + character (len=64) materialName,elemSetName + + allocate (mesh_mapFEtoCPelem(2,mesh_NcpElems)) ; mesh_mapFEtoCPelem = 0_pInt + +610 FORMAT(A300) + + cpElem = 0_pInt + rewind(unit) + do + read (unit,610,END=660) line + pos = IO_stringPos(line,maxNchunks) + select case ( IO_lc(IO_stringValue(line,pos,1)) ) + case('*material') + materialName = IO_extractValue(IO_lc(IO_stringValue(line,pos,2)),'name') ! extract name=value + materialFound = materialName /= '' ! valid name? + case('*user') + if (IO_lc(IO_stringValue(line,pos,2)) == 'material' .and. materialFound) then + do i = 1,mesh_Nmaterials ! look thru material names + if (materialName == mesh_nameMaterial(i)) exit ! found one + enddo + if (i <= mesh_Nmaterials) then ! found one? + elemSetName = mesh_mapMaterial(i) ! take corresponding elemSet + do i = 1,mesh_NelemSets ! look thru all elemSet definitions + if (elemSetName == mesh_nameElemSet(i)) then ! matched? + do j = 1,mesh_mapElemSet(1,i) + cpElem = cpElem + 1_pInt + mesh_mapFEtoCPelem(1,cpElem) = mesh_mapElemSet(1+j,i) ! store FE id + mesh_mapFEtoCPelem(2,cpElem) = cpElem ! store our id + enddo + endif + enddo + endif + materialFound = .false. + endif + endselect + enddo + +660 call qsort(mesh_mapFEtoCPelem,1,size(mesh_mapFEtoCPelem,2)) ! should be mesh_NcpElems + + if (size(mesh_mapFEtoCPelem) < 2) call IO_error(907) + + return + endsubroutine + + +!******************************************************************** +! get maximum count of nodes, IPs, IP neighbors, and subNodes +! among cpElements +! +! _maxNnodes, _maxNips, _maxNipNeighbors, _maxNsubNodes +!******************************************************************** +subroutine mesh_spectral_count_cpSizes () + + use prec, only: pInt + implicit none + + integer(pInt) t + + t = FE_mapElemtype('C3D8R') ! fake 3D hexahedral 8 node 1 IP element + + mesh_maxNnodes = FE_Nnodes(t) + mesh_maxNips = FE_Nips(t) + mesh_maxNipNeighbors = FE_NipNeighbors(t) + mesh_maxNsubNodes = FE_NsubNodes(t) + + endsubroutine + + +!******************************************************************** +! get maximum count of nodes, IPs, IP neighbors, and subNodes +! among cpElements +! +! _maxNnodes, _maxNips, _maxNipNeighbors, _maxNsubNodes +!******************************************************************** +subroutine mesh_marc_count_cpSizes (unit) + + use prec, only: pInt + use IO + implicit none + + integer(pInt), parameter :: maxNchunks = 2 + integer(pInt), dimension (1+2*maxNchunks) :: pos + character(len=300) line + + integer(pInt) unit,i,t,e + + mesh_maxNnodes = 0_pInt + mesh_maxNips = 0_pInt + mesh_maxNipNeighbors = 0_pInt + mesh_maxNsubNodes = 0_pInt + +610 FORMAT(A300) + rewind(unit) + do + read (unit,610,END=630) line + pos = IO_stringPos(line,maxNchunks) + if( IO_lc(IO_stringValue(line,pos,1)) == 'connectivity' ) then + read (unit,610,END=630) line ! Garbage line + do i=1,mesh_Nelems ! read all elements + read (unit,610,END=630) line + pos = IO_stringPos(line,maxNchunks) ! limit to id and type + e = mesh_FEasCP('elem',IO_intValue(line,pos,1)) + if (e /= 0) then + t = FE_mapElemtype(IO_stringValue(line,pos,2)) + mesh_maxNnodes = max(mesh_maxNnodes,FE_Nnodes(t)) + mesh_maxNips = max(mesh_maxNips,FE_Nips(t)) + mesh_maxNipNeighbors = max(mesh_maxNipNeighbors,FE_NipNeighbors(t)) + mesh_maxNsubNodes = max(mesh_maxNsubNodes,FE_NsubNodes(t)) + call IO_skipChunks(unit,FE_NoriginalNodes(t)-(pos(1)-2)) ! read on if FE_Nnodes exceeds node count present on current line + endif + enddo + exit + endif + enddo + +630 return + endsubroutine + + +!******************************************************************** +! get maximum count of nodes, IPs, IP neighbors, and subNodes +! among cpElements +! +! _maxNnodes, _maxNips, _maxNipNeighbors, _maxNsubNodes +!******************************************************************** + subroutine mesh_abaqus_count_cpSizes (unit) + + use prec, only: pInt + use IO + implicit none + + integer(pInt), parameter :: maxNchunks = 2 + integer(pInt), dimension (1+2*maxNchunks) :: pos + character(len=300) line + + integer(pInt) unit,i,count,t + logical inPart + + mesh_maxNnodes = 0_pInt + mesh_maxNips = 0_pInt + mesh_maxNipNeighbors = 0_pInt + mesh_maxNsubNodes = 0_pInt + +610 FORMAT(A300) + + inPart = .false. + rewind(unit) + do + read (unit,610,END=620) line + pos = IO_stringPos(line,maxNchunks) + if ( IO_lc(IO_stringValue(line,pos,1)) == '*part' ) inPart = .true. + if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. & + IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false. + + if( (inPart .or. noPart) .and. & + IO_lc(IO_stringValue(line,pos,1)) == '*element' .and. & + ( IO_lc(IO_stringValue(line,pos,2)) /= 'output' .and. & + IO_lc(IO_stringValue(line,pos,2)) /= 'matrix' .and. & + IO_lc(IO_stringValue(line,pos,2)) /= 'response' ) & + ) then + t = FE_mapElemtype(IO_extractValue(IO_lc(IO_stringValue(line,pos,2)),'type')) ! remember elem type + if (t==0) call IO_error(ID=910,ext_msg='mesh_abaqus_count_cpSizes') + count = IO_countDataLines(unit) + do i = 1,count + backspace(unit) + enddo + do i = 1,count + read (unit,610,END=620) line + pos = IO_stringPos(line,maxNchunks) ! limit to 64 nodes max + if (mesh_FEasCP('elem',IO_intValue(line,pos,1)) /= 0) then ! disregard non CP elems + mesh_maxNnodes = max(mesh_maxNnodes,FE_Nnodes(t)) + mesh_maxNips = max(mesh_maxNips,FE_Nips(t)) + mesh_maxNipNeighbors = max(mesh_maxNipNeighbors,FE_NipNeighbors(t)) + mesh_maxNsubNodes = max(mesh_maxNsubNodes,FE_NsubNodes(t)) + endif + enddo + endif + enddo + +620 return + + endsubroutine + + +!******************************************************************** +! store x,y,z coordinates of all nodes in mesh +! +! allocate globals: +! _node +!******************************************************************** + subroutine mesh_spectral_build_nodes (unit) + + use prec, only: pInt + use IO + implicit none + + integer(pInt), parameter :: maxNchunks = 7 + integer(pInt), dimension (1+2*maxNchunks) :: pos + integer(pInt) a,b,c,n,i + real(pReal) x,y,z + logical gotResolution,gotDimension + + integer(pInt) unit + character(len=64) tag + character(len=1024) line + + allocate ( mesh_node (3,mesh_Nnodes) ); mesh_node = 0_pInt + + a = 1_pInt + b = 1_pInt + c = 1_pInt + x = 1.0_pReal + y = 1.0_pReal + z = 1.0_pReal + + gotResolution = .false. + gotDimension = .false. + + rewind(unit) + do + read(unit,'(a1024)',END=100) line + if (IO_isBlank(line)) cycle ! skip empty lines + pos = IO_stringPos(line,maxNchunks) + + select case ( IO_lc(IO_StringValue(line,pos,1)) ) + case ('resolution') + gotResolution = .true. + do i = 2,6,2 + tag = IO_lc(IO_stringValue(line,pos,i)) + select case (tag) + case('a') + a = 1 + IO_intValue(line,pos,i+1) + case('b') + b = 1 + IO_intValue(line,pos,i+1) + case('c') + c = 1 + IO_intValue(line,pos,i+1) + end select + enddo + case ('dimension') + gotDimension = .true. + do i = 2,6,2 + tag = IO_lc(IO_stringValue(line,pos,i)) + select case (tag) + case('x') + x = IO_floatValue(line,pos,i+1) + case('y') + y = IO_floatValue(line,pos,i+1) + case('z') + z = IO_floatValue(line,pos,i+1) + end select + enddo + end select + if (gotDimension .and. gotResolution) exit + enddo + +! --- sanity checks --- + + if (.not. gotDimension .or. .not. gotResolution) call IO_error(42) + if (a < 2 .or. b < 2 .or. c < 2) call IO_error(43) + if (x <= 0.0_pReal .or. y <= 0.0_pReal .or. z <= 0.0_pReal) call IO_error(44) + + forall (n = 0:mesh_Nnodes-1) + mesh_node(1,n+1) = x * dble(mod(n,a) / (a-1.0_pReal)) + mesh_node(2,n+1) = y * dble(mod(n/a,b) / (b-1.0_pReal)) + mesh_node(3,n+1) = z * dble(mod(n/a/b,c) / (c-1.0_pReal)) + end forall + +100 return + + endsubroutine + + +!******************************************************************** +! store x,y,z coordinates of all nodes in mesh +! +! allocate globals: +! _node +!******************************************************************** + subroutine mesh_marc_build_nodes (unit) + + use prec, only: pInt + use IO + implicit none + + integer(pInt), dimension(5), parameter :: node_ends = (/0,10,30,50,70/) + integer(pInt), parameter :: maxNchunks = 1 + integer(pInt), dimension (1+2*maxNchunks) :: pos + character(len=300) line + + integer(pInt) unit,i,j,m + + allocate ( mesh_node (3,mesh_Nnodes) ); mesh_node = 0_pInt + +610 FORMAT(A300) + + rewind(unit) + do + read (unit,610,END=670) line + pos = IO_stringPos(line,maxNchunks) + if( IO_lc(IO_stringValue(line,pos,1)) == 'coordinates' ) then + read (unit,610,END=670) line ! skip crap line + do i=1,mesh_Nnodes + read (unit,610,END=670) line + m = mesh_FEasCP('node',IO_fixedIntValue(line,node_ends,1)) + forall (j = 1:3) mesh_node(j,m) = IO_fixedNoEFloatValue (line,node_ends,j+1) + enddo + exit + endif + enddo + +670 return + + endsubroutine + + +!******************************************************************** +! store x,y,z coordinates of all nodes in mesh +! +! allocate globals: +! _node +!******************************************************************** + subroutine mesh_abaqus_build_nodes (unit) + + use prec, only: pInt + use IO + implicit none + + integer(pInt), parameter :: maxNchunks = 4 + integer(pInt), dimension (1+2*maxNchunks) :: pos + character(len=300) line + + integer(pInt) unit,i,j,m,count + logical inPart + + allocate ( mesh_node (3,mesh_Nnodes) ); mesh_node = 0_pInt + +610 FORMAT(A300) + + inPart = .false. + rewind(unit) + do + read (unit,610,END=670) line + pos = IO_stringPos(line,maxNchunks) + if ( IO_lc(IO_stringValue(line,pos,1)) == '*part' ) inPart = .true. + if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. & + IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false. + + if( (inPart .or. noPart) .and. & + IO_lc(IO_stringValue(line,pos,1)) == '*node' .and. & + ( IO_lc(IO_stringValue(line,pos,2)) /= 'output' .and. & + IO_lc(IO_stringValue(line,pos,2)) /= 'print' .and. & + IO_lc(IO_stringValue(line,pos,2)) /= 'file' .and. & + IO_lc(IO_stringValue(line,pos,2)) /= 'response' ) & + ) then + count = IO_countDataLines(unit) ! how many nodes are defined here? + do i = 1,count + backspace(unit) ! rewind to first entry + enddo + do i = 1,count + read (unit,610,END=670) line + pos = IO_stringPos(line,maxNchunks) + m = mesh_FEasCP('node',IO_intValue(line,pos,1)) + forall (j=1:3) mesh_node(j,m) = IO_floatValue(line,pos,j+1) + enddo + endif + enddo + +670 if (size(mesh_node,2) /= mesh_Nnodes) call IO_error(909) + return + + endsubroutine + + +!******************************************************************** +! store FEid, type, mat, tex, and node list per element +! +! allocate globals: +! _element +!******************************************************************** + subroutine mesh_spectral_build_elements (unit) + + use prec, only: pInt + use IO + implicit none + + integer(pInt), parameter :: maxNchunks = 7 + integer(pInt), dimension (1+2*maxNchunks) :: pos + integer(pInt) a,b,c,e,i,homog + logical gotResolution,gotDimension,gotHomogenization + + integer(pInt) unit + character(len=1024) line + + a = 1_pInt + b = 1_pInt + c = 1_pInt + + gotResolution = .false. + gotDimension = .false. + gotHomogenization = .false. + + rewind(unit) + do + read(unit,'(a1024)',END=100) line + if (IO_isBlank(line)) cycle ! skip empty lines + pos = IO_stringPos(line,maxNchunks) + + select case ( IO_lc(IO_StringValue(line,pos,1)) ) + case ('dimension') + gotDimension = .true. + + case ('homogenization') + gotHomogenization = .true. + homog = IO_intValue(line,pos,2) + + case ('resolution') + gotResolution = .true. + do i = 2,6,2 + select case (IO_lc(IO_stringValue(line,pos,i))) + case('a') + a = IO_intValue(line,pos,i+1) + case('b') + b = IO_intValue(line,pos,i+1) + case('c') + c = IO_intValue(line,pos,i+1) + end select + enddo + end select + if (gotDimension .and. gotHomogenization .and. gotResolution) exit + enddo + +100 allocate (mesh_element (4+mesh_maxNnodes,mesh_NcpElems)) ; mesh_element = 0_pInt + + e = 0_pInt + do while (e < mesh_NcpElems) + read(unit,'(a1024)',END=110) line + if (IO_isBlank(line)) cycle ! skip empty lines + pos(1:1+2*1) = IO_stringPos(line,1) + + e = e+1 ! valid element entry + mesh_element ( 1,e) = e ! FE id + mesh_element ( 2,e) = FE_mapElemtype('C3D8R') ! elem type + mesh_element ( 3,e) = homog ! homogenization + mesh_element ( 4,e) = IO_IntValue(line,pos,1) ! microstructure + mesh_element ( 5,e) = e + (e-1)/a + (e-1)/a/b*(a+1) ! base node + mesh_element ( 6,e) = mesh_element ( 5,e) + 1 + mesh_element ( 7,e) = mesh_element ( 5,e) + (a+1) + 1 + mesh_element ( 8,e) = mesh_element ( 5,e) + (a+1) + mesh_element ( 9,e) = mesh_element ( 5,e) + (a+1)*(b+1) ! second floor base node + mesh_element (10,e) = mesh_element ( 9,e) + 1 + mesh_element (11,e) = mesh_element ( 9,e) + (a+1) + 1 + mesh_element (12,e) = mesh_element ( 9,e) + (a+1) + mesh_maxValStateVar(1) = max(mesh_maxValStateVar(1),mesh_element(3,e)) !needed for statistics + mesh_maxValStateVar(2) = max(mesh_maxValStateVar(2),mesh_element(4,e)) + enddo + +110 return + + endsubroutine + + + +!******************************************************************** +! store FEid, type, mat, tex, and node list per element +! +! allocate globals: +! _element +!******************************************************************** + subroutine mesh_marc_build_elements (unit) + + use prec, only: pInt + use IO + implicit none + + integer(pInt), parameter :: maxNchunks = 66 + integer(pInt), dimension (1+2*maxNchunks) :: pos + character(len=300) line + + integer(pInt), dimension(1+mesh_NcpElems) :: contInts + integer(pInt) unit,i,j,sv,val,e + + allocate (mesh_element (4+mesh_maxNnodes,mesh_NcpElems)) ; mesh_element = 0_pInt + +610 FORMAT(A300) + + rewind(unit) + do + read (unit,610,END=620) line + pos(1:1+2*1) = IO_stringPos(line,1) + if( IO_lc(IO_stringValue(line,pos,1)) == 'connectivity' ) then + read (unit,610,END=620) line ! Garbage line + do i = 1,mesh_Nelems + read (unit,610,END=620) line + pos = IO_stringPos(line,maxNchunks) ! limit to 64 nodes max (plus ID, type) + e = mesh_FEasCP('elem',IO_intValue(line,pos,1)) + if (e /= 0) then ! disregard non CP elems + mesh_element(1,e) = IO_IntValue (line,pos,1) ! FE id + mesh_element(2,e) = FE_mapElemtype(IO_StringValue(line,pos,2)) ! elem type + forall (j = 1:FE_Nnodes(mesh_element(2,e))) & + mesh_element(j+4,e) = IO_IntValue(line,pos,j+2) ! copy FE ids of nodes + call IO_skipChunks(unit,FE_NoriginalNodes(mesh_element(2,e))-(pos(1)-2)) ! read on if FE_Nnodes exceeds node count present on current line + endif + enddo + exit + endif + enddo + +620 rewind(unit) ! just in case "initial state" apears before "connectivity" + read (unit,610,END=620) line + do + pos(1:1+2*2) = IO_stringPos(line,2) + if( (IO_lc(IO_stringValue(line,pos,1)) == 'initial') .and. & + (IO_lc(IO_stringValue(line,pos,2)) == 'state') ) then + if (initialcondTableStyle == 2) read (unit,610,END=620) line ! read extra line for new style + read (unit,610,END=630) line ! read line with index of state var + pos(1:1+2*1) = IO_stringPos(line,1) + sv = IO_IntValue(line,pos,1) ! figure state variable index + if( (sv == 2).or.(sv == 3) ) then ! only state vars 2 and 3 of interest + read (unit,610,END=620) line ! read line with value of state var + pos(1:1+2*1) = IO_stringPos(line,1) + do while (scan(IO_stringValue(line,pos,1),'+-',back=.true.)>1) ! is noEfloat value? + val = NINT(IO_fixedNoEFloatValue(line,(/0,20/),1)) ! state var's value + mesh_maxValStateVar(sv-1) = max(val,mesh_maxValStateVar(sv-1)) ! remember max val of homogenization and microstructure index + if (initialcondTableStyle == 2) then + read (unit,610,END=630) line ! read extra line + read (unit,610,END=630) line ! read extra line + endif + contInts = IO_continousIntValues(unit,mesh_Nelems,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets) ! get affected elements + do i = 1,contInts(1) + e = mesh_FEasCP('elem',contInts(1+i)) + mesh_element(1+sv,e) = val + enddo + if (initialcondTableStyle == 0) read (unit,610,END=620) line ! ignore IP range for old table style + read (unit,610,END=630) line + pos(1:1+2*1) = IO_stringPos(line,1) + enddo + endif + else + read (unit,610,END=630) line + endif + enddo + +630 return + + endsubroutine + + + +!******************************************************************** +! store FEid, type, mat, tex, and node list per element +! +! allocate globals: +! _element +!******************************************************************** + subroutine mesh_abaqus_build_elements (unit) + + use prec, only: pInt + use IO + implicit none + + integer(pInt), parameter :: maxNchunks = 65 + integer(pInt), dimension (1+2*maxNchunks) :: pos + + integer(pInt) unit,i,j,count,e,t,homog,micro + logical inPart,materialFound + character (len=64) materialName,elemSetName + character(len=300) line + + allocate (mesh_element (4+mesh_maxNnodes,mesh_NcpElems)) ; mesh_element = 0_pInt + +610 FORMAT(A300) + + inPart = .false. + rewind(unit) + do + read (unit,610,END=620) line + pos(1:1+2*2) = IO_stringPos(line,2) + if ( IO_lc(IO_stringValue(line,pos,1)) == '*part' ) inPart = .true. + if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. & + IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false. + + if( (inPart .or. noPart) .and. & + IO_lc(IO_stringValue(line,pos,1)) == '*element' .and. & + ( IO_lc(IO_stringValue(line,pos,2)) /= 'output' .and. & + IO_lc(IO_stringValue(line,pos,2)) /= 'matrix' .and. & + IO_lc(IO_stringValue(line,pos,2)) /= 'response' ) & + ) then + t = FE_mapElemtype(IO_extractValue(IO_lc(IO_stringValue(line,pos,2)),'type')) ! remember elem type + if (t==0) call IO_error(ID=910,ext_msg='mesh_abaqus_build_elements') + count = IO_countDataLines(unit) + do i = 1,count + backspace(unit) + enddo + do i = 1,count + read (unit,610,END=620) line + pos = IO_stringPos(line,maxNchunks) ! limit to 64 nodes max + e = mesh_FEasCP('elem',IO_intValue(line,pos,1)) + if (e /= 0) then ! disregard non CP elems + mesh_element(1,e) = IO_intValue(line,pos,1) ! FE id + mesh_element(2,e) = t ! elem type + forall (j=1:FE_Nnodes(t)) & + mesh_element(4+j,e) = IO_intValue(line,pos,1+j) ! copy FE ids of nodes to position 5: + call IO_skipChunks(unit,FE_NoriginalNodes(t)-(pos(1)-1)) ! read on (even multiple lines) if FE_NoriginalNodes exceeds required node count + endif + enddo + endif + enddo + + +620 rewind(unit) ! just in case "*material" definitions apear before "*element" + + materialFound = .false. + do + read (unit,610,END=630) line + pos = IO_stringPos(line,maxNchunks) + select case ( IO_lc(IO_StringValue(line,pos,1))) + case('*material') + materialName = IO_extractValue(IO_lc(IO_StringValue(line,pos,2)),'name') ! extract name=value + materialFound = materialName /= '' ! valid name? + case('*user') + if ( IO_lc(IO_StringValue(line,pos,2)) == 'material' .and. & + materialFound ) then + do i = 1,mesh_Nmaterials ! look thru material names + if (materialName == mesh_nameMaterial(i)) exit ! found one + enddo + if (i <= mesh_Nmaterials) then ! found one? + elemSetName = mesh_mapMaterial(i) ! take corresponding elemSet + read (unit,610,END=630) line ! read homogenization and microstructure + pos(1:1+2*2) = IO_stringPos(line,2) + homog = NINT(IO_floatValue(line,pos,1)) + micro = NINT(IO_floatValue(line,pos,2)) + do i = 1,mesh_NelemSets ! look thru all elemSet definitions + if (elemSetName == mesh_nameElemSet(i)) then ! matched? + do j = 1,mesh_mapElemSet(1,i) + e = mesh_FEasCP('elem',mesh_mapElemSet(1+j,i)) + mesh_element(3,e) = homog ! store homogenization + mesh_element(4,e) = micro ! store microstructure + mesh_maxValStateVar(1) = max(mesh_maxValStateVar(1),homog) + mesh_maxValStateVar(2) = max(mesh_maxValStateVar(2),micro) + enddo + endif + enddo + endif + materialFound = .false. + endif + endselect + enddo + +630 return + + endsubroutine + + +!******************************************************************** +! get maximum count of shared elements among cpElements and +! build list of elements shared by each node in mesh +! +! _maxNsharedElems +! _sharedElem +!******************************************************************** +subroutine mesh_build_sharedElems() + +use prec, only: pInt +implicit none + +integer(pint) e, & ! element index + t, & ! element type + node, & ! CP node index + j, & ! node index per element + dim, & ! dimension index + nodeTwin ! node twin in the specified dimension +integer(pInt), dimension (mesh_Nnodes) :: node_count +integer(pInt), dimension (:), allocatable :: node_seen + +allocate(node_seen(maxval(FE_Nnodes))) + +node_count = 0_pInt + +do e = 1,mesh_NcpElems + t = mesh_element(2,e) + + node_seen = 0_pInt ! reset node duplicates + do j = 1,FE_Nnodes(t) ! check each node of element + node = mesh_FEasCP('node',mesh_element(4+j,e)) + if (all(node_seen /= node)) then + node_count(node) = node_count(node) + 1_pInt ! if FE node not yet encountered -> count it + do dim = 1,3 ! check in each dimension... + nodeTwin = mesh_nodeTwins(dim,node) + if (nodeTwin > 0) & ! if i am a twin of some node... + node_count(nodeTwin) = node_count(nodeTwin) + 1_pInt ! -> count me again for the twin node + enddo + endif + node_seen(j) = node ! remember this node to be counted already + enddo +enddo + +mesh_maxNsharedElems = maxval(node_count) ! most shared node + +allocate(mesh_sharedElem(1+mesh_maxNsharedElems,mesh_Nnodes)) +mesh_sharedElem = 0_pInt + +do e = 1,mesh_NcpElems + t = mesh_element(2,e) + node_seen = 0_pInt + do j = 1,FE_Nnodes(t) + node = mesh_FEasCP('node',mesh_element(4+j,e)) + if (all(node_seen /= node)) then + mesh_sharedElem(1,node) = mesh_sharedElem(1,node) + 1_pInt ! count for each node the connected elements + mesh_sharedElem(mesh_sharedElem(1,node)+1,node) = e ! store the respective element id + do dim = 1,3 ! check in each dimension... + nodeTwin = mesh_nodeTwins(dim,node) + if (nodeTwin > 0) then ! if i am a twin of some node... + mesh_sharedElem(1,nodeTwin) = mesh_sharedElem(1,nodeTwin) + 1_pInt ! ...count me again for the twin + mesh_sharedElem(mesh_sharedElem(1,nodeTwin)+1,nodeTwin) = e ! store the respective element id + endif + enddo + endif + node_seen(j) = node + enddo +enddo + +deallocate(node_seen) + +endsubroutine + + + +!*********************************************************** +! build up of IP neighborhood +! +! allocate globals +! _ipNeighborhood +!*********************************************************** +subroutine mesh_build_ipNeighborhood() + +use prec, only: pInt +implicit none + +integer(pInt) myElem, & ! my CP element index + myIP, & + myType, & ! my element type + myFace, & + neighbor, & ! neighor index + neighboringIPkey, & ! positive integer indicating the neighboring IP (for intra-element) and negative integer indicating the face towards neighbor (for neighboring element) + candidateIP, & + neighboringType, & ! element type of neighbor + NlinkedNodes, & ! number of linked nodes + twin_of_linkedNode, & ! node twin of a specific linkedNode + NmatchingNodes, & ! number of matching nodes + dir, & ! direction of periodicity + matchingElem, & ! CP elem number of matching element + matchingFace, & ! face ID of matching element + k, a, anchor +integer(pInt), dimension(FE_maxmaxNnodesAtIP) :: & + linkedNodes, & + matchingNodes +logical checkTwins + +allocate(mesh_ipNeighborhood(2,mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) +mesh_ipNeighborhood = 0_pInt + +linkedNodes = 0_pInt + +do myElem = 1,mesh_NcpElems ! loop over cpElems + myType = mesh_element(2,myElem) ! get elemType + do myIP = 1,FE_Nips(myType) ! loop over IPs of elem + do neighbor = 1,FE_NipNeighbors(myType) ! loop over neighbors of IP + neighboringIPkey = FE_ipNeighbor(neighbor,myIP,myType) + + !*** if the key is positive, the neighbor is inside the element + !*** that means, we have already found our neighboring IP + + if (neighboringIPkey > 0) then + mesh_ipNeighborhood(1,neighbor,myIP,myElem) = myElem + mesh_ipNeighborhood(2,neighbor,myIP,myElem) = neighboringIPkey + + + !*** if the key is negative, the neighbor resides in a neighboring element + !*** that means, we have to look through the face indicated by the key and see which element is behind that face + + elseif (neighboringIPkey < 0) then ! neighboring element's IP + myFace = -neighboringIPkey + call mesh_faceMatch(myElem, myFace, matchingElem, matchingFace) ! get face and CP elem id of face match + if (matchingElem > 0_pInt) then ! found match? + neighboringType = mesh_element(2,matchingElem) + + !*** trivial solution if neighbor has only one IP + + if (FE_Nips(neighboringType) == 1_pInt) then + mesh_ipNeighborhood(1,neighbor,myIP,myElem) = matchingElem + mesh_ipNeighborhood(2,neighbor,myIP,myElem) = 1_pInt + cycle + endif + + !*** find those nodes which build the link to the neighbor + + NlinkedNodes = 0_pInt + linkedNodes = 0_pInt + do a = 1,FE_maxNnodesAtIP(myType) ! figure my anchor nodes on connecting face + anchor = FE_nodesAtIP(a,myIP,myType) + if (anchor /= 0_pInt) then ! valid anchor node + if (any(FE_nodeOnFace(:,myFace,myType) == anchor)) then ! ip anchor sits on face? + NlinkedNodes = NlinkedNodes + 1_pInt + linkedNodes(NlinkedNodes) = mesh_element(4+anchor,myElem) ! FE id of anchor node + else ! something went wrong with the linkage, since not all anchors sit on my face + NlinkedNodes = 0_pInt + linkedNodes = 0_pInt + exit + endif + endif + enddo + + !*** loop through the ips of my neighbor + !*** and try to find an ip with matching nodes + !*** also try to match with node twins + +checkCandidateIP: do candidateIP = 1,FE_Nips(neighboringType) + NmatchingNodes = 0_pInt + matchingNodes = 0_pInt + do a = 1,FE_maxNnodesAtIP(neighboringType) ! check each anchor node of that ip + anchor = FE_nodesAtIP(a,candidateIP,neighboringType) + if (anchor /= 0_pInt) then ! valid anchor node + if (any(FE_nodeOnFace(:,matchingFace,neighboringType) == anchor)) then ! sits on matching face? + NmatchingNodes = NmatchingNodes + 1_pInt + matchingNodes(NmatchingNodes) = mesh_element(4+anchor,matchingElem) ! FE id of neighbor's anchor node + else ! no matching, because not all nodes sit on the matching face + NmatchingNodes = 0_pInt + matchingNodes = 0_pInt + exit + endif + endif + enddo + + if (NmatchingNodes /= NlinkedNodes) & ! this ip has wrong count of anchors on face + cycle checkCandidateIP + + !*** check "normal" nodes whether they match or not + + checkTwins = .false. + do a = 1,NlinkedNodes + if (all(matchingNodes /= linkedNodes(a))) then ! this linkedNode does not match any matchingNode + checkTwins = .true. + exit ! no need to search further + endif + enddo + + !*** if no match found, then also check node twins + + if(checkTwins) then +periodicityDirection: do dir = 1,3 + do a = 1,NlinkedNodes + twin_of_linkedNode = mesh_nodeTwins(dir,linkedNodes(a)) + if (twin_of_linkedNode == 0_pInt & ! twin of linkedNode does not exist... + .or. all(matchingNodes /= twin_of_linkedNode)) then ! ... or it does not match any matchingNode + if (dir < 3) then ! no match in this direction... + cycle periodicityDirection ! ... so try in different direction + else ! no matching in any direction... + cycle checkCandidateIP ! ... so check next candidateIP + endif + endif + enddo + exit periodicityDirection + enddo periodicityDirection + endif + + !*** we found a match !!! + + mesh_ipNeighborhood(1,neighbor,myIP,myElem) = matchingElem + mesh_ipNeighborhood(2,neighbor,myIP,myElem) = candidateIP + exit checkCandidateIP + + enddo checkCandidateIP + endif ! end of valid external matching + endif ! end of internal/external matching + enddo + enddo +enddo + +endsubroutine + + + +!*********************************************************** +! assignment of coordinates for subnodes in each cp element +! +! allocate globals +! _subNodeCoord +!*********************************************************** + subroutine mesh_build_subNodeCoords() + + use prec, only: pInt,pReal + implicit none + + integer(pInt) e,t,n,p + + allocate(mesh_subNodeCoord(3,mesh_maxNnodes+mesh_maxNsubNodes,mesh_NcpElems)) ; mesh_subNodeCoord = 0.0_pReal + + do e = 1,mesh_NcpElems ! loop over cpElems + t = mesh_element(2,e) ! get elemType + do n = 1,FE_Nnodes(t) + mesh_subNodeCoord(:,n,e) = mesh_node(:,mesh_FEasCP('node',mesh_element(4+n,e))) ! loop over nodes of this element type + enddo + do n = 1,FE_NsubNodes(t) ! now for the true subnodes + do p = 1,FE_Nips(t) ! loop through possible parent nodes + if (FE_subNodeParent(p,n,t) > 0) & ! valid parent node + mesh_subNodeCoord(:,n+FE_Nnodes(t),e) = & + mesh_subNodeCoord(:,n+FE_Nnodes(t),e) + & + mesh_node(:,mesh_FEasCP('node',mesh_element(4+FE_subNodeParent(p,n,t),e))) ! add up parents + enddo + mesh_subNodeCoord(:,n+FE_Nnodes(t),e) = mesh_subNodeCoord(:,n+FE_Nnodes(t),e) / count(FE_subNodeParent(:,n,t) > 0) + enddo + enddo + + return + + endsubroutine + + +!*********************************************************** +! calculation of IP volume +! +! allocate globals +! _ipVolume +!*********************************************************** + subroutine mesh_build_ipVolumes() + + use prec, only: pInt + use math, only: math_volTetrahedron + implicit none + + integer(pInt) e,f,t,i,j,k,n + integer(pInt), parameter :: Ntriangles = FE_NipFaceNodes-2 ! each interface is made up of this many triangles + logical(pInt), dimension(mesh_maxNnodes+mesh_maxNsubNodes) :: gravityNode ! flagList to find subnodes determining center of grav + real(pReal), dimension(3,mesh_maxNnodes+mesh_maxNsubNodes) :: gravityNodePos ! coordinates of subnodes determining center of grav + real(pReal), dimension(3,FE_NipFaceNodes) :: nPos ! coordinates of nodes on IP face + real(pReal), dimension(Ntriangles,FE_NipFaceNodes) :: volume ! volumes of possible tetrahedra + real(pReal), dimension(3) :: centerOfGravity + + allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems)) ; mesh_ipVolume = 0.0_pReal + allocate(mesh_ipCenterOfGravity(3,mesh_maxNips,mesh_NcpElems)) ; mesh_ipCenterOfGravity = 0.0_pReal + + do e = 1,mesh_NcpElems ! loop over cpElems + t = mesh_element(2,e) ! get elemType + do i = 1,FE_Nips(t) ! loop over IPs of elem + gravityNode = .false. ! reset flagList + gravityNodePos = 0.0_pReal ! reset coordinates + do f = 1,FE_NipNeighbors(t) ! loop over interfaces of IP + do n = 1,FE_NipFaceNodes ! loop over nodes on interface + gravityNode(FE_subNodeOnIPFace(n,f,i,t)) = .true. + gravityNodePos(:,FE_subNodeOnIPFace(n,f,i,t)) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e) + enddo + enddo + + do j = 1,mesh_maxNnodes+mesh_maxNsubNodes-1 ! walk through entire flagList except last + if (gravityNode(j)) then ! valid node index + do k = j+1,mesh_maxNnodes+mesh_maxNsubNodes ! walk through remainder of list + if (gravityNode(k) .and. all(abs(gravityNodePos(:,j) - gravityNodePos(:,k)) < 1.0e-36_pReal)) then ! found duplicate + gravityNode(j) = .false. ! delete first instance + gravityNodePos(:,j) = 0.0_pReal + exit ! continue with next suspect + endif + enddo + endif + enddo + centerOfGravity = sum(gravityNodePos,2)/count(gravityNode) + + do f = 1,FE_NipNeighbors(t) ! loop over interfaces of IP and add tetrahedra which connect to CoG + forall (n = 1:FE_NipFaceNodes) nPos(:,n) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e) + forall (n = 1:FE_NipFaceNodes, j = 1:Ntriangles) & ! start at each interface node and build valid triangles to cover interface + volume(j,n) = math_volTetrahedron(nPos(:,n), & ! calc volume of respective tetrahedron to CoG + nPos(:,1+mod(n+j-1,FE_NipFaceNodes)), & + nPos(:,1+mod(n+j-0,FE_NipFaceNodes)), & + centerOfGravity) + mesh_ipVolume(i,e) = mesh_ipVolume(i,e) + sum(volume) ! add contribution from this interface + enddo + mesh_ipVolume(i,e) = mesh_ipVolume(i,e) / FE_NipFaceNodes ! renormalize with interfaceNodeNum due to loop over them + mesh_ipCenterOfGravity(:,i,e) = centerOfGravity + enddo + enddo + return + + endsubroutine + + +!*********************************************************** +! calculation of IP interface areas +! +! allocate globals +! _ipArea, _ipAreaNormal +!*********************************************************** + subroutine mesh_build_ipAreas() + + use prec, only: pInt,pReal + use math + implicit none + + integer(pInt) e,f,t,i,j,n + integer(pInt), parameter :: Ntriangles = FE_NipFaceNodes-2 ! each interface is made up of this many triangles + real(pReal), dimension (3,FE_NipFaceNodes) :: nPos ! coordinates of nodes on IP face + real(pReal), dimension(3,Ntriangles,FE_NipFaceNodes) :: normal + real(pReal), dimension(Ntriangles,FE_NipFaceNodes) :: area + + allocate(mesh_ipArea(mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) ; mesh_ipArea = 0.0_pReal + allocate(mesh_ipAreaNormal(3,mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) ; mesh_ipAreaNormal = 0.0_pReal + + do e = 1,mesh_NcpElems ! loop over cpElems + t = mesh_element(2,e) ! get elemType + do i = 1,FE_Nips(t) ! loop over IPs of elem + do f = 1,FE_NipNeighbors(t) ! loop over interfaces of IP + forall (n = 1:FE_NipFaceNodes) nPos(:,n) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e) + forall (n = 1:FE_NipFaceNodes, j = 1:Ntriangles) ! start at each interface node and build valid triangles to cover interface + normal(:,j,n) = math_vectorproduct(nPos(:,1+mod(n+j-1,FE_NipFaceNodes)) - nPos(:,n), & ! calc their normal vectors + nPos(:,1+mod(n+j-0,FE_NipFaceNodes)) - nPos(:,n)) + area(j,n) = sqrt(sum(normal(:,j,n)*normal(:,j,n))) ! and area + end forall + forall (n = 1:FE_NipFaceNodes, j = 1:Ntriangles, area(j,n) > 0.0_pReal) & + normal(:,j,n) = normal(:,j,n) / area(j,n) ! make unit normal + + mesh_ipArea(f,i,e) = sum(area) / (FE_NipFaceNodes*2.0_pReal) ! area of parallelograms instead of triangles + mesh_ipAreaNormal(:,f,i,e) = sum(sum(normal,3),2) / count(area > 0.0_pReal) ! average of all valid normals + enddo + enddo + enddo + return + + endsubroutine + + +!*********************************************************** +! assignment of twin nodes for each cp node +! +! allocate globals +! _nodeTwins +!*********************************************************** +subroutine mesh_build_nodeTwins() + +use prec, only: pInt, pReal +implicit none + +integer(pInt) dir, & ! direction of periodicity + node, & + minimumNode, & + maximumNode, & + n1, & + n2 +integer(pInt), dimension(mesh_Nnodes+1) :: minimumNodes, maximumNodes ! list of surface nodes (minimum and maximum coordinate value) with first entry giving the number of nodes +real(pReal) minCoord, maxCoord, & ! extreme positions in one dimension + tolerance ! tolerance below which positions are assumed identical +real(pReal), dimension(3) :: distance ! distance between two nodes in all three coordinates +logical, dimension(mesh_Nnodes) :: node_seen + +allocate(mesh_nodeTwins(3,mesh_Nnodes)) +mesh_nodeTwins = 0_pInt + +tolerance = 0.001_pReal * minval(mesh_ipVolume) ** 0.333_pReal + +do dir = 1,3 ! check periodicity in directions of x,y,z + if (mesh_periodicSurface(dir)) then ! only if periodicity is requested + + + !*** find out which nodes sit on the surface + !*** and have a minimum or maximum position in this dimension + + minimumNodes = 0_pInt + maximumNodes = 0_pInt + minCoord = minval(mesh_node(dir,:)) + maxCoord = maxval(mesh_node(dir,:)) + do node = 1,mesh_Nnodes ! loop through all nodes and find surface nodes + if (abs(mesh_node(dir,node) - minCoord) <= tolerance) then + minimumNodes(1) = minimumNodes(1) + 1_pInt + minimumNodes(minimumNodes(1)+1) = node + elseif (abs(mesh_node(dir,node) - maxCoord) <= tolerance) then + maximumNodes(1) = maximumNodes(1) + 1_pInt + maximumNodes(maximumNodes(1)+1) = node + endif + enddo + + + !*** find the corresponding node on the other side with the same position in this dimension + + node_seen = .false. + do n1 = 1,minimumNodes(1) + minimumNode = minimumNodes(n1+1) + if (node_seen(minimumNode)) & + cycle + do n2 = 1,maximumNodes(1) + maximumNode = maximumNodes(n2+1) + distance = abs(mesh_node(:,minimumNode) - mesh_node(:,maximumNode)) + if (sum(distance) - distance(dir) <= tolerance) then ! minimum possible distance (within tolerance) + mesh_nodeTwins(dir,minimumNode) = maximumNode + mesh_nodeTwins(dir,maximumNode) = minimumNode + node_seen(maximumNode) = .true. ! remember this node, we don't have to look for his partner again + exit + endif + enddo + enddo + + endif +enddo + +endsubroutine + + + +!*********************************************************** +! write statistics regarding input file parsing +! to the output file +! +!*********************************************************** +subroutine mesh_tell_statistics() + +use prec, only: pInt +use math, only: math_range +use IO, only: IO_error +use debug, only: verboseDebugger + +implicit none + +integer(pInt), dimension (:,:), allocatable :: mesh_HomogMicro +character(len=64) fmt + +integer(pInt) i,e,n,f,t + +if (mesh_maxValStateVar(1) < 1_pInt) call IO_error(110) ! no homogenization specified +if (mesh_maxValStateVar(2) < 1_pInt) call IO_error(120) ! no microstructure specified + +allocate (mesh_HomogMicro(mesh_maxValStateVar(1),mesh_maxValStateVar(2))); mesh_HomogMicro = 0_pInt +do e = 1,mesh_NcpElems + if (mesh_element(3,e) < 1_pInt) call IO_error(110,e) ! no homogenization specified + if (mesh_element(4,e) < 1_pInt) call IO_error(120,e) ! no homogenization specified + mesh_HomogMicro(mesh_element(3,e),mesh_element(4,e)) = & + mesh_HomogMicro(mesh_element(3,e),mesh_element(4,e)) + 1 ! count combinations of homogenization and microstructure +enddo + +if (verboseDebugger) then + !$OMP CRITICAL (write2out) + write (6,*) + write (6,*) "Input Parser: SUBNODE COORDINATES" + write (6,*) + write(6,'(a5,x,a5,x,a15,x,a15,x,a20,3(x,a12))') 'elem','IP','IP neighbor','IPFaceNodes','subNodeOnIPFace','x','y','z' + do e = 1,mesh_NcpElems ! loop over cpElems + t = mesh_element(2,e) ! get elemType + do i = 1,FE_Nips(t) ! loop over IPs of elem + do f = 1,FE_NipNeighbors(t) ! loop over interfaces of IP + do n = 1,FE_NipFaceNodes ! loop over nodes on interface + write(6,'(i5,x,i5,x,i15,x,i15,x,i20,3(x,f12.8))') e,i,f,n,FE_subNodeOnIPFace(n,f,i,t),& + mesh_subNodeCoord(1,FE_subNodeOnIPFace(n,f,i,t),e),& + mesh_subNodeCoord(2,FE_subNodeOnIPFace(n,f,i,t),e),& + mesh_subNodeCoord(3,FE_subNodeOnIPFace(n,f,i,t),e) + enddo + enddo + enddo + enddo + write(6,*) + write(6,*) 'Input Parser: IP COORDINATES' + write(6,'(a5,x,a5,3(x,a12))') 'elem','IP','x','y','z' + do e = 1,mesh_NcpElems + do i = 1,FE_Nips(mesh_element(2,e)) + write (6,'(i5,x,i5,3(x,f12.8))') e, i, mesh_ipCenterOfGravity(:,i,e) + enddo + enddo + write (6,*) + write (6,*) "Input Parser: ELEMENT VOLUME" + write (6,*) + write (6,"(a13,x,e15.8)") "total volume", sum(mesh_ipVolume) + write (6,*) + write (6,"(a5,x,a5,x,a15,x,a5,x,a15,x,a16)") "elem","IP","volume","face","area","-- normal --" + do e = 1,mesh_NcpElems + do i = 1,FE_Nips(mesh_element(2,e)) + write (6,"(i5,x,i5,x,e15.8)") e,i,mesh_IPvolume(i,e) + do f = 1,FE_NipNeighbors(mesh_element(2,e)) + write (6,"(i33,x,e15.8,x,3(f6.3,x))") f,mesh_ipArea(f,i,e),mesh_ipAreaNormal(:,f,i,e) + enddo + enddo + enddo + write (6,*) + write (6,*) "Input Parser: NODE TWINS" + write (6,*) + write(6,'(a6,3(3(x),a6))') ' node','twin_x','twin_y','twin_z' + do n = 1,mesh_Nnodes ! loop over cpNodes + write(6,'(i6,3(3(x),i6))') n, mesh_nodeTwins(1:3,n) + enddo + write(6,*) + write(6,*) "Input Parser: IP NEIGHBORHOOD" + write(6,*) + write(6,"(a10,x,a10,x,a10,x,a3,x,a13,x,a13)") "elem","IP","neighbor","","elemNeighbor","ipNeighbor" + do e = 1,mesh_NcpElems ! loop over cpElems + t = mesh_element(2,e) ! get elemType + do i = 1,FE_Nips(t) ! loop over IPs of elem + do n = 1,FE_NipNeighbors(t) ! loop over neighbors of IP + write (6,"(i10,x,i10,x,i10,x,a3,x,i13,x,i13)") e,i,n,'-->',mesh_ipNeighborhood(1,n,i,e),mesh_ipNeighborhood(2,n,i,e) + enddo + enddo + enddo + !$OMP END CRITICAL (write2out) +endif + + +!$OMP CRITICAL (write2out) + write (6,*) + write (6,*) "Input Parser: STATISTICS" + write (6,*) + write (6,*) mesh_Nelems, " : total number of elements in mesh" + write (6,*) mesh_NcpElems, " : total number of CP elements in mesh" + write (6,*) mesh_Nnodes, " : total number of nodes in mesh" + write (6,*) mesh_maxNnodes, " : max number of nodes in any CP element" + write (6,*) mesh_maxNips, " : max number of IPs in any CP element" + write (6,*) mesh_maxNipNeighbors, " : max number of IP neighbors in any CP element" + write (6,*) mesh_maxNsubNodes, " : max number of (additional) subnodes in any CP element" + write (6,*) mesh_maxNsharedElems, " : max number of CP elements sharing a node" + write (6,*) + write (6,*) "Input Parser: HOMOGENIZATION/MICROSTRUCTURE" + write (6,*) + write (6,*) mesh_maxValStateVar(1), " : maximum homogenization index" + write (6,*) mesh_maxValStateVar(2), " : maximum microstructure index" + write (6,*) + write (fmt,"(a,i5,a)") "(9(x),a2,x,",mesh_maxValStateVar(2),"(i8))" + write (6,fmt) "+-",math_range(mesh_maxValStateVar(2)) + write (fmt,"(a,i5,a)") "(i8,x,a2,x,",mesh_maxValStateVar(2),"(i8))" + do i=1,mesh_maxValStateVar(1) ! loop over all (possibly assigned) homogenizations + write (6,fmt) i,"| ",mesh_HomogMicro(i,:) ! loop over all (possibly assigned) microstrcutures + enddo + write(6,*) + write(6,*) "Input Parser: ADDITIONAL MPIE OPTIONS" + write(6,*) + write(6,*) "periodic surface : ", mesh_periodicSurface + write(6,*) + call flush(6) +!$OMP END CRITICAL (write2out) + +deallocate(mesh_HomogMicro) + +endsubroutine + + +END MODULE mesh + diff --git a/code/mpie_spectral_single.f90 b/code/mpie_spectral_single.f90 index e4eaa5ea8..8f0ee34a8 100644 --- a/code/mpie_spectral_single.f90 +++ b/code/mpie_spectral_single.f90 @@ -561,7 +561,7 @@ program mpie_spectral enddo ! end looping over steps in current loadcase enddo ! end looping over loadcases close(538) - call dfftw_destroy_plan(plan_fft(1)); call dfftw_destroy_plan(plan_fft(2)) + call sfftw_destroy_plan(plan_fft(1)); call sfftw_destroy_plan(plan_fft(2)) end program mpie_spectral diff --git a/processing/post/spectral_post_single.py b/processing/post/spectral_post_single.py new file mode 100644 index 000000000..4bab00158 --- /dev/null +++ b/processing/post/spectral_post_single.py @@ -0,0 +1,368 @@ +#!/usr/bin/python +# -*- coding: iso-8859-1 -*- + +# This script is used for the post processing of the results achieved by the spectral method. +# As it reads in the data coming from "materialpoint_results", it can be adopted to the data +# computed using the FEM solvers. Until now, its capable to handle elements with one IP in a regular order +# written by M. Diehl, m.diehl@mpie.de + +import os,sys,re,array,struct,numpy, time, postprocessingMath, math + +class vector: + x,y,z = [None,None,None] + + def __init__(self,coords): + self.x = coords[0] + self.y = coords[1] + self.z = coords[2] + +class element: + items = [] + type = None + + def __init__(self,nodes,type): + self.items = nodes + self.type = type + +class element_scalar: + id = None + value = None + + def __init__(self,node,value): + self.id = node + self.value = value + + +class MPIEspectral_result: + + file = None + dataOffset = 0 + N_elemental_scalars = 0 + resolution = [0,0,0] + dimension = [0.0,0.0,0.0] + theTitle = '' + wd = '' + extrapolate = '' + N_increments = 0 + increment = 0 + N_nodes = 0 + N_node_scalars = 0 + N_elements = 0 + N_element_scalars = 0 + N_element_tensors = 0 + theNodes = [] + theElements = [] + + def __init__(self,filename): + + self.file = open(filename, 'rb') + + self.title = self._keyedString('load') + self.wd = self._keyedString('workingdir') + self.geometry = self._keyedString('geometry') + self.N_increments = self._keyedInt('increments') + self.N_element_scalars = self._keyedInt('materialpoint_sizeResults') + self.resolution = self._keyedPackedArray('resolution',3,'i') + self.N_nodes = (self.resolution[0]+1)*(self.resolution[1]+1)*(self.resolution[2]+1) + self.N_elements = self.resolution[0]*self.resolution[1]*self.resolution[2] + self.dimension = self._keyedPackedArray('dimension',3,'f') + a = self.resolution[0]+1 + b = self.resolution[1]+1 + c = self.resolution[2]+1 + self.file.seek(0) + self.dataOffset = self.file.read(2048).find('eoh')+7 + + + def __str__(self): + return '\n'.join([ + 'title: %s'%self.title, + 'workdir: %s'%self.wd, + 'extrapolation: %s'%self.extrapolate, + 'increments: %i'%self.N_increments, + 'increment: %i'%self.increment, + 'nodes: %i'%self.N_nodes, + 'resolution: %s'%(','.join(map(str,self.resolution))), + 'dimension: %s'%(','.join(map(str,self.dimension))), + 'elements: %i'%self.N_elements, + 'nodal_scalars: %i'%self.N_node_scalars, + 'elemental scalars: %i'%self.N_element_scalars, + 'end of header: %i'%self.dataOffset, + ] + ) + + def _keyedPackedArray(self,identifier,length = 3,type = 'f'): + match = {'f': 4,'i': 4} #correct??? + self.file.seek(0) + m = re.search('%s%s'%(identifier,'(.{%i})'%(match[type])*length),self.file.read(2048),re.DOTALL) + values = [] + if m: + for i in m.groups(): + values.append(struct.unpack(type,i)[0]) + return values + + def _keyedInt(self,identifier): + value = None + self.file.seek(0) + m = re.search('%s%s'%(identifier,'(.{4})'),self.file.read(2048),re.DOTALL) + if m: + value = struct.unpack('i',m.group(1))[0] + return value + + def _keyedString(self,identifier): + value = None + self.file.seek(0) + m = re.search(r'(.{4})%s(.*?)\1'%identifier,self.file.read(2048),re.DOTALL) + if m: + value = m.group(2) + return value + + def extrapolation(self,value): + self.extrapolate = value + + def element_scalar(self,elem,idx): + self.file.seek(self.dataOffset+(self.increment*(4+self.N_elements*self.N_element_scalars*4+4) + 4+(elem*self.N_element_scalars + idx)*4)) + value = struct.unpack('f',self.file.read(4))[0] + return [elemental_scalar(node,value) for node in self.theElements[elem].items] + +def readScalar(resolution,file,distance,startingPosition,offset): + currentPosition = startingPosition+offset*4+4 - distance*4 # we add distance later on + field = numpy.zeros([resolution[0],resolution[1],resolution[2]], 'f') + for z in range(0,resolution[2]): + for y in range(0,resolution[1]): + for x in range(0,resolution[0]): + currentPosition = currentPosition + distance*4 + p.file.seek(currentPosition) + field[x][y][z]=struct.unpack('f',p.file.read(4))[0] + return field + +def readTensor(resolution,file,distance,startingPosition,offset): + currentPosition = startingPosition+offset*4+4 - distance*4 # we add distance later on + field = numpy.zeros([resolution[0],resolution[1],resolution[2],3,3], 'f') + for z in range(0,resolution[2]): + for y in range(0,resolution[1]): + for x in range(0,resolution[0]): + currentPosition = currentPosition + distance*4 + p.file.seek(currentPosition) + for i in range(0,3): + for j in range(0,3): + field[x][y][z][i][j]=struct.unpack('f',p.file.read(4))[0] + return field + +#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +def calculateCauchyStress(p_stress,defgrad,res): +#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + c_stress = numpy.zeros([res[0],res[1],res[2],3,3],'f') + for z in range(res[2]): + for y in range(res[1]): + for x in range(res[0]): + jacobi = numpy.linalg.det(defgrad[x,y,z]) + c_stress[x,y,z] = numpy.dot(p_stress[x,y,z],numpy.transpose(defgrad[x,y,z]))/jacobi + return c_stress + + +# function writes scalar values to a mesh (geometry) +def writeVtkAscii(filename,geometry,scalar,resolution): + prodnn=(p.resolution[0]+1)*(p.resolution[1]+1)*(p.resolution[2]+1) + vtk = open(filename, 'w') + vtk.write('# vtk DataFile Version 3.1\n') # header + vtk.write('just a test\n') # header + vtk.write('ASCII\n') # header + vtk.write('DATASET UNSTRUCTURED_GRID\n') # header + vtk.write('POINTS ') # header + vtk.write(str(prodnn)) # header + vtk.write(' FLOAT\n') # header + +# nodes + for k in range (resolution[2]+1): + for j in range (resolution[1]+1): + for i in range (resolution[0]+1): + vtk.write('\t'.join(map(str,geometry[i,j,k]))+'\n') + vtk.write('\n') + vtk.write('CELLS ') + vtk.write(str(resolution[0]*resolution[1]*resolution[2])) + vtk.write('\t') + vtk.write(str(resolution[0]*resolution[1]*resolution[2]*9)) + vtk.write('\n') + +# elements + for i in range (resolution[2]): + for j in range (resolution[1]): + for k in range (resolution[0]): + vtk.write('8') + vtk.write('\t') + base = i*(resolution[1]+1)*(resolution[2]+1)+j*(resolution[1]+1)+k + vtk.write(str(base)) + vtk.write('\t') + vtk.write(str(base+1)) + vtk.write('\t') + vtk.write(str(base+resolution[1]+2)) + vtk.write('\t') + vtk.write(str(base+resolution[1]+1)) + vtk.write('\t') + base = base + (resolution[1]+1)*(resolution[2]+1) + vtk.write(str(base)) + vtk.write('\t') + vtk.write(str(base+1)) + vtk.write('\t') + vtk.write(str(base+resolution[1]+2)) + vtk.write('\t') + vtk.write(str(base+resolution[1]+1)) + vtk.write('\n') + vtk.write('\n') + vtk.write('CELL_TYPES ') + vtk.write('\t') + vtk.write(str(resolution[0]*resolution[1]*resolution[2])) + vtk.write('\n') + for i in range (resolution[0]*resolution[1]*resolution[2]): + vtk.write('12\n') + vtk.write('\nCELL_DATA ') # header + vtk.write(str(resolution[0]*resolution[1]*resolution[2])) # header + vtk.write('\n') # header + vtk.write('SCALARS HorizontalSpeed float\n') # header + vtk.write('LOOKUP_TABLE default\n') # header + for k in range (resolution[2]): + for j in range (resolution[1]): + for i in range (resolution[0]): + vtk.write(str(scalar[i,j,k])) + vtk.write('\n') + return + +# function writes scalar values to a point field +def writeVtkAsciiDots(filename,coordinates,scalar,resolution): + prodnn=(p.resolution[0])*(p.resolution[1])*(p.resolution[2]) + vtk = open(filename, 'w') + vtk.write('# vtk DataFile Version 3.1\n') # header + vtk.write('just a test\n') # header + vtk.write('ASCII\n') # header + vtk.write('DATASET UNSTRUCTURED_GRID\n') # header + vtk.write('POINTS ') # header + vtk.write(str(prodnn)) # header + vtk.write(' FLOAT\n') # header +# points + for k in range (resolution[2]): + for j in range (resolution[1]): + for i in range (resolution[0]): + vtk.write('\t'.join(map(str,coordinates[i,j,k]))+'\n') + vtk.write('\n') + vtk.write('CELLS ') + vtk.write(str(prodnn)) + vtk.write('\t') + vtk.write(str(prodnn*2)) + vtk.write('\n') + for i in range(prodnn): + vtk.write('1\t' + str(i) + '\n') + vtk.write('CELL_TYPES ') + vtk.write('\t') + vtk.write(str(prodnn)) + vtk.write('\n') + for i in range (prodnn): + vtk.write('1\n') + vtk.write('\nPOINT_DATA ') # header + vtk.write(str(prodnn)) # header + vtk.write('\n') # header + vtk.write('SCALARS HorizontalSpeed float\n') # header + vtk.write('LOOKUP_TABLE default\n') # header + for k in range (resolution[2]): + for j in range (resolution[1]): + for i in range (resolution[0]): + vtk.write(str(scalar[i,j,k])) + vtk.write('\n') + return + +# functiongives the corner box for the average defgrad +def writeVtkAsciidefgrad_av(filename,diag,defgrad): + + points = numpy.array([\ + [0.0,0.0,0.0,],\ + [diag[0],0.0,0.0,],\ + [diag[0],diag[1],0.0,],\ + [0.0,diag[1],0.0,],\ + [0.0,0.0,diag[2],],\ + [diag[0],0.0,diag[2],],\ + [diag[0],diag[1],diag[2],],\ + [0.0,diag[1],diag[2],]]\ + ) + vtk = open(filename, 'w') + vtk.write('# vtk DataFile Version 3.1\n') # header + vtk.write('just a test\n') # header + vtk.write('ASCII\n') # header + vtk.write('DATASET UNSTRUCTURED_GRID\n') # header + vtk.write('POINTS 8') # header + vtk.write(' FLOAT\n') # header + + # points + for p in range (8): + vtk.write('\t'.join(map(str,numpy.dot(defgrad_av,points[p])))+'\n') + vtk.write('\n') + vtk.write('CELLS 8 16\n') + vtk.write('\n'.join(['1\t%i'%i for i in range(8)])+'\n') + vtk.write('CELL_TYPES 8\n') + vtk.write('\n'.join(['1']*8)+'\n') + + return + +print '*********************************************************************************' +print 'Post Processing for Material subroutine for BVP solution using spectral method' +print '*********************************************************************************\n' + +#reading in the header of the results file +name = 'dipl32_shear' +p = MPIEspectral_result(name+'.spectralOut') +p.extrapolation('') +print p + +# Ended reading of header + +res_x=p.resolution[0] +res_y=p.resolution[1] +res_z=p.resolution[2] + + +ms=numpy.zeros([res_x,res_y,res_z,3], 'f') +print 'data structure' +for i in range(p.N_element_scalars): + c_pos = p.dataOffset + i*4.0 + 4.0 + p.file.seek(c_pos) + print(i, struct.unpack('f',p.file.read(4))) +for i in range(1,2): + c_pos = p.dataOffset + i*(p.N_element_scalars*4*p.N_elements + 8) + defgrad = readTensor(p.resolution,p.file,p.N_element_scalars,c_pos,7) + defgrad_av = postprocessingMath.tensor_avg(res_x,res_y,res_z,defgrad) + centroids_coord = postprocessingMath.deformed(res_x,res_y,res_z,p.dimension,defgrad,defgrad_av) + undeformed = postprocessingMath.mesh(p.resolution[0],p.resolution[1],p.resolution[2],p.dimension,defgrad_av,centroids_coord) + #writeVtkAscii(name+'-mesh-undeformed.vtk',undeformed,defgrad[:,:,:,1,2],p.resolution) + +for i in range(240,241): + c_pos = p.dataOffset + i*(p.N_element_scalars*4*p.N_elements + 8) #8 accounts for header&footer + defgrad = readTensor(p.resolution,p.file,p.N_element_scalars,c_pos,7) + defgrad_av = postprocessingMath.tensor_avg(res_x,res_y,res_z,defgrad) + #defgrad = numpy.zeros([p.resolution[0],p.resolution[1],p.resolution[2],3,3], 'f') + #for z in range(p.resolution[2]): + # defgrad[:,:,z,1,2] = (2.0*z)/(p.resolution[2]-1.0)+ 3.8*math.sin(z*20.0/(p.resolution[2]-1.0)*2*math.pi) + # defgrad[:,:,z,0,2] = (2.0*z)/(p.resolution[2]-1.0)+ 5.0*math.cos(z/(p.resolution[2]-1.0)*2*math.pi) + #defgrad[:,:,:,0,0] = 1.0 + #defgrad[:,:,:,1,1] = 1.0 + #defgrad[:,:,:,2,2] = 1.0 + #logstrain = postprocessingMath.logstrain_mat(res_x,res_y,res_z,defgrad) + #logstrain2 = postprocessingMath.logstrain_spat(res_x,res_y,res_z,defgrad) + #p_stress = readTensor(p.resolution,p.file,p.N_element_scalars,c_pos,52) + #c_stress = postprocessingMath.calculate_cauchy(res_x,res_y,res_z,defgrad,p_stress) + #vm = postprocessingMath.calculate_mises(res_x,res_y,res_z,c_stress) + #defgrad_av = postprocessingMath.tensor_avg(res_x,res_y,res_z,defgrad) + #subroutine inverse_reconstruction(res_x,res_y,res_z,reference_configuration,current_configuration,defgrad) + centroids_coord = postprocessingMath.deformed(res_x,res_y,res_z,p.dimension,defgrad,defgrad_av) + deformed = postprocessingMath.mesh(p.resolution[0],p.resolution[1],p.resolution[2],p.dimension,defgrad_av,centroids_coord) + #writeVtkAscii(name+'-mesh-deformed.vtk',deformed,defgrad[:,:,:,1,2],p.resolution) + defgrad = postprocessingMath.inverse_reconstruction(res_x,res_y,res_z,undeformed,deformed) + defgrad_av = postprocessingMath.tensor_avg(res_x,res_y,res_z,defgrad) + centroids_coord = postprocessingMath.deformed(res_x,res_y,res_z,p.dimension,defgrad,defgrad_av) + centroids_coord2 = postprocessingMath.deformed_fft(res_x,res_y,res_z,p.dimension,defgrad,defgrad_av,1.0) + ms = postprocessingMath.mesh(p.resolution[0],p.resolution[1],p.resolution[2],p.dimension,defgrad_av,centroids_coord) + ms2 = postprocessingMath.mesh(p.resolution[0],p.resolution[1],p.resolution[2],p.dimension,defgrad_av,centroids_coord2) + writeVtkAscii(name+'-mesh-usual-%s.vtk'%i,ms,defgrad[:,:,:,1,2],p.resolution) + writeVtkAscii(name+'-mesh-fft-%s.vtk'%i,ms2,defgrad[:,:,:,1,2],p.resolution) + #writeVtkAsciidefgrad_av(name+'-box-%i.vtk'%i,p.dimension,defgrad_av) + #writeVtkAsciiDots(name+'-points-%i.vtk'%i,centroids_coord,grain,p.resolution) + sys.stdout.flush() + +