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

@@ -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

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
+      wallFunc = 4.0_pReal * sqrt(2.0_pReal) / 3.0_pReal * sqrt(1.0_pReal - tauRel) / tauRel
-      boltzmannProbability = dexp(- constitutive_nonlocal_Qeff0(s,myInstance) * wallFunc / (kB * Temperature))
+      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

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

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)))
+         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)/dsqrt(math_mul3x3(lattice_fcc_systemSlip(4:6,i),lattice_fcc_systemSlip(4:6,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)))
+         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)/dsqrt(math_mul3x3(lattice_fcc_systemTwin(4:6,i),lattice_fcc_systemTwin(4:6,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)))
+         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)/dsqrt(math_mul3x3(lattice_bcc_systemSlip(4:6,i),lattice_bcc_systemSlip(4:6,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)))
+         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)/dsqrt(math_mul3x3(lattice_bcc_systemTwin(4:6,i),lattice_bcc_systemTwin(4:6,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))
+         sd(:,i) = hex_d/sqrt(math_mul3x3(hex_d,hex_d))
-         sn(:,i) = hex_n/dsqrt(math_mul3x3(hex_n,hex_n))
+         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))
+         td(:,i) = hex_d/sqrt(math_mul3x3(hex_d,hex_d))
-         tn(:,i) = hex_n/dsqrt(math_mul3x3(hex_n,hex_n))
+         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

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
+precision :=double
-	ifort -openmp -o cpspectral.exe mpie_spectral.o CPFEM.a libfftw3_threads.a libfftw3.a constitutive.a advanced.a basics.a -lpthread
+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

code/makefile_single

@@ -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

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(1) = cos(2.0_pReal*pi*rnd(1))*sqrt(rnd(3))
- math_qRnd(2) = dsin(2.0_pReal*pi*rnd(2))*dsqrt(1.0_pReal-rnd(3))
+ math_qRnd(2) = sin(2.0_pReal*pi*rnd(2))*sqrt(1.0_pReal-rnd(3))
- math_qRnd(3) = dcos(2.0_pReal*pi*rnd(2))*dsqrt(1.0_pReal-rnd(3))
+ math_qRnd(3) = cos(2.0_pReal*pi*rnd(2))*sqrt(1.0_pReal-rnd(3))
- math_qRnd(4) = dsin(2.0_pReal*pi*rnd(1))*dsqrt(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(1) = 0.5_pReal * sqrt(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(2) = 0.5_pReal * sqrt(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(3) = 0.5_pReal * sqrt(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(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))
+ C1 = cos(Euler(1))
- C = dcos(Euler(2))
+ C = cos(Euler(2))
- C2 = dcos(Euler(3))
+ C2 = cos(Euler(3))
- S1 = dsin(Euler(1))
+ S1 = sin(Euler(1))
- S = dsin(Euler(2))
+ S = sin(Euler(2))
- S2 = dsin(Euler(3))
+ 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))
+ c = cos(halfangles(2))
- s = dsin(halfangles(2))
+ s = sin(halfangles(2))
  
- math_EulerToQuaternion(1) = dcos(halfangles(1)+halfangles(3)) * c
+ math_EulerToQuaternion(1) = cos(halfangles(1)+halfangles(3)) * c
- math_EulerToQuaternion(2) = dcos(halfangles(1)-halfangles(3)) * s
+ math_EulerToQuaternion(2) = cos(halfangles(1)-halfangles(3)) * s
- math_EulerToQuaternion(3) = dsin(halfangles(1)-halfangles(3)) * s
+ math_EulerToQuaternion(3) = sin(halfangles(1)-halfangles(3)) * s
- math_EulerToQuaternion(4) = dsin(halfangles(1)+halfangles(3)) * c
+ 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)
+   s = sin(omega)
-   c = dcos(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)
+   s = sin(omega/2.0_pReal)
-   c = dcos(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
+ if (abs(math_QuaternionToEuler(2)) < 1.0e-3_pReal) then
-   math_QuaternionToEuler(1) = 2.0_pReal*dacos(Q(1))
+   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
+ halfAngle = acos(max(-1.0_pReal, min(1.0_pReal, Q(1))))            ! limit to [-1,1] --> 0 to 180 deg
- sinHalfAngle = dsin(halfAngle)
+ 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)
+   PHI=acos(arg)
-   Y1=2*RHO**(1.0_pReal/3.0_pReal)*DCOS(PHI/3.0_pReal)
+   Y1=2*RHO**(1.0_pReal/3.0_pReal)*cos(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)
+   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)*DCOS(PHI/3.0_pReal+4.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

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

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
 

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()
+
+