diff --git a/.gitignore b/.gitignore index c34f2f0b7..3f1fc3458 100644 --- a/.gitignore +++ b/.gitignore @@ -1,5 +1,6 @@ *.pyc *.hdf5 +*.xdmf *.bak *~ bin diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index 0c3a6248a..3ee9bc1c5 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -906,7 +906,7 @@ class DADF5(): n_nodes = 8 for i in f['/geometry/T_c']: - vtk_geom.InsertNextCell(n_nodes,vtk_type,i-1) + vtk_geom.InsertNextCell(vtk_type,n_nodes,i-1) elif mode == 'Point': Points = vtk.vtkPoints() diff --git a/python/damask/orientation.py b/python/damask/orientation.py index dc6a74504..abdaba661 100644 --- a/python/damask/orientation.py +++ b/python/damask/orientation.py @@ -1078,12 +1078,12 @@ class Orientation: if isinstance(lattice, Lattice): self.lattice = lattice else: - self.lattice = Lattice(lattice) # assume string + self.lattice = Lattice(lattice) # assume string if isinstance(rotation, Rotation): self.rotation = rotation else: - self.rotation = Rotation(rotation) # assume quaternion + self.rotation = Rotation.fromQuaternion(rotation) # assume quaternion def disorientation(self, other, diff --git a/python/damask/table.py b/python/damask/table.py index 62a2540f6..582089a54 100644 --- a/python/damask/table.py +++ b/python/damask/table.py @@ -23,11 +23,7 @@ class Table(): """ self.comments = [] if comments is None else [c for c in comments] - if hasattr(data,'columns'): - self.data = data - self.data.columns = [''] * len(self.data.columns) - else: - self.data = pd.DataFrame(data=data) + self.data = pd.DataFrame(data=data) self.shapes = shapes self.__label_condensed() diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index b324a5afc..5c8266b2a 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -65,12 +65,11 @@ program DAMASK_spectral currentLoadcase = 0, & !< current load case inc, & !< current increment in current load case totalIncsCounter = 0, & !< total # of increments - fileUnit = 0, & !< file unit for reading load case and writing results - myStat, & statUnit = 0, & !< file unit for statistics output stagIter, & nActiveFields = 0 character(len=6) :: loadcase_string + character(len=pStringLen), dimension(:), allocatable :: fileContent character(len=1024) :: & incInfo type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases @@ -141,17 +140,14 @@ program DAMASK_spectral !-------------------------------------------------------------------------------------------------- ! reading information from load case file and to sanity checks + fileContent = IO_read_ASCII(trim(loadCaseFile)) + allocate (loadCases(0)) ! array of load cases - open(newunit=fileunit,iostat=myStat,file=trim(loadCaseFile),action='read') - if (myStat /= 0) call IO_error(100,el=myStat,ext_msg=trim(loadCaseFile)) - do - read(fileUnit, '(A)', iostat=myStat) line - if ( myStat /= 0) exit - if (IO_isBlank(line)) cycle ! skip empty lines - - currentLoadCase = currentLoadCase + 1 - + do currentLoadCase = 1, size(fileContent) + line = fileContent(currentLoadCase) + if (IO_isBlank(line)) cycle chunkPos = IO_stringPos(line) + do i = 1, chunkPos(1) ! reading compulsory parameters for loadcase select case (IO_lc(IO_stringValue(line,chunkPos,i))) case('l','fdot','dotf','f') @@ -300,7 +296,7 @@ program DAMASK_spectral endif reportAndCheck loadCases = [loadCases,newLoadCase] ! load case is ok, append it enddo - close(fileUnit) + !-------------------------------------------------------------------------------------------------- ! doing initialization depending on active solvers diff --git a/src/mesh/DAMASK_FEM.f90 b/src/mesh/DAMASK_FEM.f90 index 9b9b95b91..269296e84 100644 --- a/src/mesh/DAMASK_FEM.f90 +++ b/src/mesh/DAMASK_FEM.f90 @@ -66,7 +66,7 @@ program DAMASK_FEM PetscInt :: faceSet, currentFaceSet PetscInt :: field, dimPlex PetscErrorCode :: ierr - + integer(kind(COMPONENT_UNDEFINED_ID)) :: ID external :: & quit @@ -166,45 +166,28 @@ program DAMASK_FEM !-------------------------------------------------------------------------------------------------- ! boundary condition information - case('x') ! X displacement field - do field = 1, nActiveFields - if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_MECH_ID) then - do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents - if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_MECH_X_ID) then - loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = & - .true. - loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = & - IO_floatValue(line,chunkPos,i+1) - endif - enddo - endif - enddo - case('y') ! Y displacement field - do field = 1, nActiveFields - if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_MECH_ID) then - do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents - if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_MECH_Y_ID) then - loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = & - .true. - loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = & - IO_floatValue(line,chunkPos,i+1) - endif - enddo - endif - enddo - case('z') ! Z displacement field - do field = 1, nActiveFields - if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_MECH_ID) then - do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents - if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_MECH_Z_ID) then - loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = & - .true. - loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = & - IO_floatValue(line,chunkPos,i+1) - endif - enddo - endif - enddo + case('x','y','z') + select case(IO_lc(IO_stringValue(line,chunkPos,i))) + case('x') + ID = COMPONENT_MECH_X_ID + case('y') + ID = COMPONENT_MECH_Y_ID + case('z') + ID = COMPONENT_MECH_Z_ID + end select + + do field = 1, nActiveFields + if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_MECH_ID) then + do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents + if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == ID) then + loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = & + .true. + loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = & + IO_floatValue(line,chunkPos,i+1) + endif + enddo + endif + enddo end select enddo; enddo close(fileUnit) diff --git a/src/mesh/FEM_mech.f90 b/src/mesh/FEM_mech.f90 index 3a19f67c7..0280629ea 100644 --- a/src/mesh/FEM_mech.f90 +++ b/src/mesh/FEM_mech.f90 @@ -67,28 +67,33 @@ contains subroutine FEM_mech_init(fieldBC) type(tFieldBC), intent(in) :: fieldBC + DM :: mech_mesh PetscFE :: mechFE PetscQuadrature :: mechQuad, functional PetscDS :: mechDS PetscDualSpace :: mechDualSpace + DMLabel, dimension(:),pointer :: pBCLabel DMLabel :: BCLabel + PetscInt, dimension(:), pointer :: pNumComp, pNumDof, pBcField, pBcPoint - PetscInt :: numBC, bcSize, nc + PetscInt :: numBC, bcSize, nc, & + field, faceSet, topologDim, nNodalPoints, & + cellStart, cellEnd, cell, basis + IS :: bcPoint - IS, pointer :: pBcComps(:), pBcPoints(:) + IS, dimension(:), pointer :: pBcComps, pBcPoints PetscSection :: section - PetscInt :: field, faceSet, topologDim, nNodalPoints + PetscReal, dimension(:), pointer :: qPointsP, qWeightsP, & - nodalPointsP, nodalWeightsP - PetscReal, allocatable, target :: nodalPoints(:), nodalWeights(:) - PetscScalar, pointer :: px_scal(:) - PetscScalar, allocatable, target :: x_scal(:) + nodalPointsP, nodalWeightsP,pV0, pCellJ, pInvcellJ PetscReal :: detJ PetscReal, allocatable, target :: cellJMat(:,:) - PetscReal, pointer :: pV0(:), pCellJ(:), pInvcellJ(:) - PetscInt :: cellStart, cellEnd, cell, basis - character(len=7), parameter :: prefix = 'mechFE_' + + PetscScalar, pointer :: px_scal(:) + PetscScalar, allocatable, target :: x_scal(:) + + character(len=*), parameter :: prefix = 'mechFE_' PetscErrorCode :: ierr write(6,'(/,a)') ' <<<+- FEM_mech init -+>>>' @@ -125,13 +130,19 @@ subroutine FEM_mech_init(fieldBC) ! Setup FEM mech boundary conditions call DMGetLabel(mech_mesh,'Face Sets',BCLabel,ierr); CHKERRQ(ierr) call DMPlexLabelComplete(mech_mesh,BCLabel,ierr); CHKERRQ(ierr) +#if (PETSC_VERSION_MINOR < 11) call DMGetSection(mech_mesh,section,ierr); CHKERRQ(ierr) +#else + call DMGetLocalSection(mech_mesh,section,ierr); CHKERRQ(ierr) +#endif allocate(pnumComp(1), source=dimPlex) - allocate(pnumDof(dimPlex+1), source = 0) + allocate(pnumDof(0:dimPlex), source = 0) do topologDim = 0, dimPlex call DMPlexGetDepthStratum(mech_mesh,topologDim,cellStart,cellEnd,ierr) CHKERRQ(ierr) - call PetscSectionGetDof(section,cellStart,pnumDof(topologDim+1),ierr) + call PetscSectionGetDof(section,cellStart,pnumDof(topologDim),ierr) + write(6,*) 'start',cellStart,'end',cellEnd + write(6,*) 'topologDim',topologDim,'numDOF',pNumDOF(topologDim) CHKERRQ(ierr) enddo numBC = 0 @@ -163,9 +174,15 @@ subroutine FEM_mech_init(fieldBC) endif endif enddo; enddo +#if (PETSC_VERSION_MINOR < 11) call DMPlexCreateSection(mech_mesh,dimPlex,1,pNumComp,pNumDof, & - numBC,pBcField,pBcComps,pBcPoints,PETSC_NULL_IS, & - section,ierr) + numBC,pBcField,pBcComps,pBcPoints,PETSC_NULL_IS,section,ierr) +#else + allocate(pBClabel(1),source=BClabel) + call DMPlexCreateSection(mech_mesh,pBClabel,pNumComp,pNumDof, & + numBC,pBcField,pBcComps,pBcPoints,PETSC_NULL_IS,section,ierr) + +#endif CHKERRQ(ierr) call DMSetSection(mech_mesh,section,ierr); CHKERRQ(ierr) do faceSet = 1, numBC @@ -196,13 +213,11 @@ subroutine FEM_mech_init(fieldBC) call VecSet(solution ,0.0,ierr); CHKERRQ(ierr) call VecSet(solution_rate ,0.0,ierr); CHKERRQ(ierr) allocate(x_scal(cellDof)) - allocate(nodalPoints (dimPlex)) - allocate(nodalWeights(1)) - nodalPointsP => nodalPoints - nodalWeightsP => nodalWeights + allocate(nodalWeightsP(1)) + allocate(nodalPointsP(dimPlex)) allocate(pv0(dimPlex)) - allocate(pcellJ(dimPlex*dimPlex)) - allocate(pinvcellJ(dimPlex*dimPlex)) + allocate(pcellJ(dimPlex**2)) + allocate(pinvcellJ(dimPlex**2)) allocate(cellJMat(dimPlex,dimPlex)) call DMGetSection(mech_mesh,section,ierr); CHKERRQ(ierr) call DMGetDS(mech_mesh,mechDS,ierr); CHKERRQ(ierr) @@ -212,7 +227,7 @@ subroutine FEM_mech_init(fieldBC) call DMPlexGetHeightStratum(mech_mesh,0,cellStart,cellEnd,ierr) CHKERRQ(ierr) do cell = cellStart, cellEnd-1 !< loop over all elements - x_scal = 0.0 + x_scal = 0.0_pReal call DMPlexComputeCellGeometryAffineFEM(mech_mesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr) CHKERRQ(ierr) cellJMat = reshape(pCellJ,shape=[dimPlex,dimPlex]) @@ -221,7 +236,7 @@ subroutine FEM_mech_init(fieldBC) CHKERRQ(ierr) call PetscQuadratureGetData(functional,dimPlex,nc,nNodalPoints,nodalPointsP,nodalWeightsP,ierr) CHKERRQ(ierr) - x_scal(basis+1:basis+dimPlex) = pV0 + matmul(transpose(cellJMat),nodalPointsP + 1.0) + x_scal(basis+1:basis+dimPlex) = pV0 + matmul(transpose(cellJMat),nodalPointsP + 1.0_pReal) enddo px_scal => x_scal call DMPlexVecSetClosure(mech_mesh,section,solution_local,cell,px_scal,INSERT_ALL_VALUES,ierr) @@ -283,6 +298,9 @@ end function FEM_mech_solution subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr) DM :: dm_local + PetscObject,intent(in) :: dummy + PetscErrorCode :: ierr + PetscDS :: prob Vec :: x_local, f_local, xx_local PetscSection :: section @@ -294,10 +312,10 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr) qPt, basis, comp, cidx PetscReal :: detFAvg PetscReal :: BMat(dimPlex*dimPlex,cellDof) - PetscObject,intent(in) :: dummy + PetscInt :: bcSize IS :: bcPoints - PetscErrorCode :: ierr + allocate(pV0(dimPlex)) allocate(pcellJ(dimPlex**2)) @@ -316,7 +334,7 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr) call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,ierr) CHKERRQ(ierr) call utilities_projectBCValues(x_local,section,0,field-1,bcPoints, & - 0.0,params%fieldBC%componentBC(field)%Value(face),params%timeinc) + 0.0_pReal,params%fieldBC%componentBC(field)%Value(face),params%timeinc) call ISDestroy(bcPoints,ierr); CHKERRQ(ierr) endif endif @@ -403,29 +421,35 @@ end subroutine FEM_mech_formResidual !-------------------------------------------------------------------------------------------------- subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr) - DM :: dm_local + + DM :: dm_local + Mat :: Jac_pre, Jac + PetscObject, intent(in) :: dummy + PetscErrorCode :: ierr + PetscDS :: prob Vec :: x_local, xx_local - Mat :: Jac_pre, Jac + PetscSection :: section, gSection + + PetscReal, dimension(1, cellDof) :: MatB + PetscReal, dimension(dimPlex**2,cellDof) :: BMat, BMatAvg, MatA + PetscReal, dimension(3,3) :: F, FAvg, FInv PetscReal :: detJ PetscReal, dimension(:), pointer :: basisField, basisFieldDer, & pV0, pCellJ, pInvcellJ + + PetscScalar, dimension(:), pointer :: pK_e, x_scal + + PetscScalar,dimension(cellDOF,cellDOF), target :: K_e + PetscScalar,dimension(cellDOF,cellDOF) :: K_eA , & + K_eB + PetscInt :: cellStart, cellEnd, cell, field, face, & qPt, basis, comp, cidx,bcSize - PetscScalar,dimension(cellDOF,cellDOF), target :: K_e, & - K_eA , & - K_eB - PetscScalar, target :: K_eVec(cellDof*cellDof) - PetscReal :: BMat (dimPlex*dimPlex,cellDof), & - BMatAvg(dimPlex*dimPlex,cellDof), & - MatA (dimPlex*dimPlex,cellDof), & - MatB (1 ,cellDof) - PetscScalar, dimension(:), pointer :: pK_e, x_scal - PetscReal, dimension(3,3) :: F, FAvg, FInv - PetscObject, intent(in) :: dummy + IS :: bcPoints - PetscErrorCode :: ierr + allocate(pV0(dimPlex)) allocate(pcellJ(dimPlex**2)) @@ -440,7 +464,7 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr) call DMGetGlobalSection(dm_local,gSection,ierr); CHKERRQ(ierr) call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr) - call VecWAXPY(x_local,1.0,xx_local,solution_local,ierr); CHKERRQ(ierr) + call VecWAXPY(x_local,1.0_pReal,xx_local,solution_local,ierr); CHKERRQ(ierr) do field = 1, dimPlex; do face = 1, mesh_Nboundaries if (params%fieldBC%componentBC(field)%Mask(face)) then call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,ierr) @@ -448,7 +472,7 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr) call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,ierr) CHKERRQ(ierr) call utilities_projectBCValues(x_local,section,0,field-1,bcPoints, & - 0.0,params%fieldBC%componentBC(field)%Value(face),params%timeinc) + 0.0_pReal,params%fieldBC%componentBC(field)%Value(face),params%timeinc) call ISDestroy(bcPoints,ierr); CHKERRQ(ierr) endif endif @@ -501,13 +525,16 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr) (matmul(matmul(transpose(BMatAvg), & reshape(FInv(1:dimPlex,1:dimPlex),shape=[dimPlex*dimPlex,1],order=[2,1])),MatB) + & K_eB)/real(dimPlex) - else K_e = K_eA endif - K_e = K_e + eps*math_identity2nd(cellDof) - K_eVec = reshape(K_e, [cellDof*cellDof])*abs(detJ) - pK_e => K_eVec + K_e = (K_e + eps*math_identity2nd(cellDof)) * abs(detJ) +#ifndef __INTEL_COMPILER + pK_e(1:cellDOF**2) => K_e +#else + ! https://software.intel.com/en-us/forums/intel-fortran-compiler/topic/782230 (bug) + allocate(pK_e(cellDOF**2),source = reshape(K_e,[cellDOF**2])) +#endif call DMPlexMatSetClosure(dm_local,section,gSection,Jac,cell,pK_e,ADD_VALUES,ierr) CHKERRQ(ierr) call DMPlexVecRestoreClosure(dm_local,section,x_local,cell,x_scal,ierr) @@ -536,18 +563,18 @@ subroutine FEM_mech_forward(guess,timeinc,timeinc_old,fieldBC) type(tFieldBC), intent(in) :: & fieldBC - real(pReal), intent(in) :: & + real(pReal), intent(in) :: & timeinc_old, & timeinc - logical, intent(in) :: & + logical, intent(in) :: & guess - PetscInt :: field, face - DM :: dm_local - Vec :: x_local - PetscSection :: section - PetscInt :: bcSize - IS :: bcPoints - PetscErrorCode :: ierr + + PetscInt :: field, face, bcSize + DM :: dm_local + Vec :: x_local + PetscSection :: section + IS :: bcPoints + PetscErrorCode :: ierr !-------------------------------------------------------------------------------------------------- ! forward last inc @@ -557,7 +584,7 @@ subroutine FEM_mech_forward(guess,timeinc,timeinc_old,fieldBC) call SNESGetDM(mech_snes,dm_local,ierr); CHKERRQ(ierr) !< retrieve mesh info from mech_snes into dm_local call DMGetSection(dm_local,section,ierr); CHKERRQ(ierr) call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr) - call VecSet(x_local,0.0,ierr); CHKERRQ(ierr) + call VecSet(x_local,0.0_pReal,ierr); CHKERRQ(ierr) call DMGlobalToLocalBegin(dm_local,solution,INSERT_VALUES,x_local,ierr) !< retrieve my partition of global solution vector CHKERRQ(ierr) call DMGlobalToLocalEnd(dm_local,solution,INSERT_VALUES,x_local,ierr) @@ -570,7 +597,7 @@ subroutine FEM_mech_forward(guess,timeinc,timeinc_old,fieldBC) call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,ierr) CHKERRQ(ierr) call utilities_projectBCValues(solution_local,section,0,field-1,bcPoints, & - 0.0,fieldBC%componentBC(field)%Value(face),timeinc_old) + 0.0_pReal,fieldBC%componentBC(field)%Value(face),timeinc_old) call ISDestroy(bcPoints,ierr); CHKERRQ(ierr) endif endif diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index 1303f0df8..cb30f1590 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -23,7 +23,7 @@ module FEM_utilities private !-------------------------------------------------------------------------------------------------- - logical, public :: cutBack = .false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill + logical, public :: cutBack = .false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill integer, public, parameter :: maxFields = 6 integer, public :: nActiveFields = 0 @@ -56,15 +56,15 @@ module FEM_utilities !-------------------------------------------------------------------------------------------------- ! derived types type, public :: tSolutionState !< return type of solution from FEM solver variants - logical :: converged = .true. - logical :: stagConverged = .true. - integer :: iterationsNeeded = 0 + logical :: converged = .true. + logical :: stagConverged = .true. + integer :: iterationsNeeded = 0 end type tSolutionState type, public :: tComponentBC integer(kind(COMPONENT_UNDEFINED_ID)) :: ID - real(pReal), allocatable :: Value(:) - logical, allocatable :: Mask(:) + real(pReal), allocatable, dimension(:) :: Value + logical, allocatable, dimension(:) :: Mask end type tComponentBC type, public :: tFieldBC @@ -79,8 +79,8 @@ module FEM_utilities outputfrequency = 1, & !< frequency of result writes logscale = 0 !< linear/logarithmic time inc flag logical :: followFormerTrajectory = .true. !< follow trajectory of former loadcase - integer, allocatable :: faceID(:) - type(tFieldBC), allocatable :: fieldBC(:) + integer, allocatable, dimension(:) :: faceID + type(tFieldBC), allocatable, dimension(:) :: fieldBC end type tLoadCase public :: & @@ -88,6 +88,7 @@ module FEM_utilities utilities_constitutiveResponse, & utilities_projectBCValues, & FIELD_MECH_ID, & + COMPONENT_UNDEFINED_ID, & COMPONENT_MECH_X_ID, & COMPONENT_MECH_Y_ID, & COMPONENT_MECH_Z_ID diff --git a/src/mesh/FEM_zoo.f90 b/src/mesh/FEM_zoo.f90 index 5ef37e462..5ab1c3102 100644 --- a/src/mesh/FEM_zoo.f90 +++ b/src/mesh/FEM_zoo.f90 @@ -37,7 +37,7 @@ contains !-------------------------------------------------------------------------------------------------- subroutine FEM_Zoo_init - write(6,'(/,a)') ' <<<+- FEM_Zoo init -+>>>' + write(6,'(/,a)') ' <<<+- FEM_Zoo init -+>>>'; flush(6) !-------------------------------------------------------------------------------------------------- ! 2D linear @@ -53,7 +53,7 @@ subroutine FEM_Zoo_init ! 2D quadratic FEM_Zoo_nQuadrature(2,2) = 3 - allocate(FEM_Zoo_QuadratureWeights(2,2)%p(3)) + allocate(FEM_Zoo_QuadratureWeights(2,2)%p(3)) FEM_Zoo_QuadratureWeights(2,2)%p(1:3) = 1.0_pReal/3.0_pReal allocate(FEM_Zoo_QuadraturePoints (2,2)%p(6)) diff --git a/src/mesh_FEM.f90 b/src/mesh_FEM.f90 index afbdaf7f9..014c2d84c 100644 --- a/src/mesh_FEM.f90 +++ b/src/mesh_FEM.f90 @@ -26,9 +26,11 @@ module mesh integer, public, protected :: & mesh_Nboundaries, & - mesh_NcpElems, & !< total number of CP elements in mesh mesh_NcpElemsGlobal + integer :: & + mesh_NcpElems !< total number of CP elements in mesh + !!!! BEGIN DEPRECATED !!!!! integer, public, protected :: & mesh_maxNips !< max number of IPs in any CP element @@ -44,7 +46,7 @@ module mesh mesh_ipVolume, & !< volume associated with IP (initially!) mesh_node0 !< node x,y,z coordinates (initially!) - real(pReal), dimension(:,:,:), allocatable, public :: & + real(pReal), dimension(:,:,:), allocatable :: & mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!) DM, public :: geomMesh @@ -177,7 +179,7 @@ subroutine mesh_init call IO_error(602,ext_msg='IP') ! selected element does not have requested IP FEsolving_execElem = [ 1,mesh_NcpElems ] ! parallel loop bounds set to comprise all DAMASK elements - FEsolving_execIP = spread([1,FE_Nips(FE_geomtype(mesh_element(2,1))),2,nElems) + FEsolving_execIP = spread([1,FE_Nips(FE_geomtype(mesh_element(2,1)))],2,mesh_NcpElems) allocate(mesh_node0(3,mesh_Nnodes),source=0.0_pReal) @@ -193,21 +195,17 @@ end subroutine mesh_init !-------------------------------------------------------------------------------------------------- subroutine mesh_FEM_build_ipVolumes(dimPlex) - PetscInt :: dimPlex + PetscInt,intent(in):: dimPlex PetscReal :: vol - PetscReal, target :: cent(dimPlex), norm(dimPlex) - PetscReal, pointer :: pCent(:), pNorm(:) + PetscReal, pointer,dimension(:) :: pCent, pNorm PetscInt :: cellStart, cellEnd, cell PetscErrorCode :: ierr - if (.not. allocated(mesh_ipVolume)) then - allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems)) - mesh_ipVolume = 0.0_pReal - endif + allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems),source=0.0_pReal) call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr) - pCent => cent - pNorm => norm + allocate(pCent(dimPlex)) + allocate(pNorm(dimPlex)) do cell = cellStart, cellEnd-1 call DMPlexComputeCellGeometryFVM(geomMesh,cell,vol,pCent,pNorm,ierr) CHKERRQ(ierr) @@ -225,8 +223,7 @@ subroutine mesh_FEM_build_ipCoordinates(dimPlex,qPoints) PetscInt, intent(in) :: dimPlex PetscReal, intent(in) :: qPoints(mesh_maxNips*dimPlex) - PetscReal, target :: v0(dimPlex), cellJ(dimPlex*dimPlex), invcellJ(dimPlex*dimPlex) - PetscReal, pointer :: pV0(:), pCellJ(:), pInvcellJ(:) + PetscReal, pointer,dimension(:) :: pV0, pCellJ, pInvcellJ PetscReal :: detJ PetscInt :: cellStart, cellEnd, cell, qPt, dirI, dirJ, qOffset PetscErrorCode :: ierr @@ -234,9 +231,9 @@ subroutine mesh_FEM_build_ipCoordinates(dimPlex,qPoints) allocate(mesh_ipCoordinates(3,mesh_maxNips,mesh_NcpElems),source=0.0_pReal) - pV0 => v0 - pCellJ => cellJ - pInvcellJ => invcellJ + allocate(pV0(dimPlex)) + allocatE(pCellJ(dimPlex**2)) + allocatE(pinvCellJ(dimPlex**2)) call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr) do cell = cellStart, cellEnd-1 !< loop over all elements call DMPlexComputeCellGeometryAffineFEM(geomMesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr)