diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index ad22c7998..df5a10f60 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -75,7 +75,7 @@ variables: MSC: "$MSC2019" IntelMarc: "$IntelCompiler17_8" IntelAbaqus: "$IntelCompiler16_4" - HDF5Marc: "HDF5/1.10.4/Intel-17.8" + HDF5Marc: "HDF5/1.10.5/Intel-17.8" # ++++++++++++ Documentation ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Doxygen1_8_15: "Documentation/Doxygen/1.8.15" # ------------ Defaults ---------------------------------------------- diff --git a/PRIVATE b/PRIVATE index cdf792813..214c69be8 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit cdf79281330b2cf6e11426685b1d0c999efc24cf +Subproject commit 214c69be8b51adb39eb7ad25b139727c8b98afce diff --git a/VERSION b/VERSION index 6edf451c2..040bd7cf0 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-860-g0dd14a06 +v2.0.3-922-gcebbfc90 diff --git a/installation/mods_MarcMentat/2019/Marc_tools/include_linux64 b/installation/mods_MarcMentat/2019/Marc_tools/include_linux64 index 0ef7c5733..290055ed3 100644 --- a/installation/mods_MarcMentat/2019/Marc_tools/include_linux64 +++ b/installation/mods_MarcMentat/2019/Marc_tools/include_linux64 @@ -63,6 +63,7 @@ else INTEGER_PATH=/$MARC_INTEGER_SIZE fi +FCOMP=ifort INTELPATH="/opt/intel/compilers_and_libraries_2017/linux" # find the root directory of the compiler installation: @@ -103,9 +104,6 @@ if test "$DAMASK_HDF5" = "ON";then H5FC="$(h5fc -shlib -show)" HDF5_LIB=${H5FC//ifort/} FCOMP="$H5FC -DDAMASK_HDF5" - echo $FCOMP -else - FCOMP=ifort fi # AEM @@ -531,7 +529,7 @@ else FORT_OPT=" $FORT_OPT -save -zero" fi if test "$MARCHDF" = "HDF"; then - FORT_OPT="$FORT_OPT -DMARCHDF=$MARCHDF $HDF_INCLUDE" + FORT_OPT="$FORT_OPT -DMARCHDF=$MARCHDF" fi FORTLOW="$FCOMP $FORT_OPT $PROFILE -O0 $I8FFLAGS -I$MARC_SOURCE/common \ @@ -757,7 +755,7 @@ SECLIBS="-L$MARC_LIB -llapi" SOLVERLIBS="${BCSSOLVERLIBS} ${VKISOLVERLIBS} ${CASISOLVERLIBS} ${MF2SOLVERLIBS} \ $MKLLIB -L$MARC_MKL -liomp5 \ - $MARC_LIB/blas_src.a ${ACSI_LIB}/ACSI_MarcLib.a $KDTREE2_LIB/kdtree2.a $HDF_LIBS $HDF_LIB" + $MARC_LIB/blas_src.a ${ACSI_LIB}/ACSI_MarcLib.a $KDTREE2_LIB/kdtree2.a $HDF5_LIB" SOLVERLIBS_DLL=${SOLVERLIBS} if test "$AEM_DLL" -eq 1 diff --git a/processing/post/DADF5_vtk_cells.py b/processing/post/DADF5_vtk_cells.py index 0b8fd07de..700ce4404 100755 --- a/processing/post/DADF5_vtk_cells.py +++ b/processing/post/DADF5_vtk_cells.py @@ -3,6 +3,7 @@ import os import argparse +import h5py import numpy as np import vtk from vtk.util import numpy_support @@ -40,21 +41,30 @@ if options.con is None: options.con=[] for filename in options.filenames: results = damask.DADF5(filename) - if results.structured: # for grid solvers use rectilinear grid - rGrid = vtk.vtkRectilinearGrid() + if results.structured: # for grid solvers use rectilinear grid + grid = vtk.vtkRectilineagrid() coordArray = [vtk.vtkDoubleArray(), vtk.vtkDoubleArray(), vtk.vtkDoubleArray(), ] - rGrid.SetDimensions(*(results.grid+1)) + grid.SetDimensions(*(results.grid+1)) for dim in [0,1,2]: for c in np.linspace(0,results.size[dim],1+results.grid[dim]): coordArray[dim].InsertNextValue(c) - rGrid.SetXCoordinates(coordArray[0]) - rGrid.SetYCoordinates(coordArray[1]) - rGrid.SetZCoordinates(coordArray[2]) + grid.SetXCoordinates(coordArray[0]) + grid.SetYCoordinates(coordArray[1]) + grid.SetZCoordinates(coordArray[2]) + else: + nodes = vtk.vtkPoints() + with h5py.File(filename) as f: + nodes.SetData(numpy_support.numpy_to_vtk(f['/geometry/x_n'][()],deep=True)) + grid = vtk.vtkUnstructuredGrid() + grid.SetPoints(nodes) + grid.Allocate(f['/geometry/T_c'].shape[0]) + for i in f['/geometry/T_c']: + grid.InsertNextCell(vtk.VTK_HEXAHEDRON,8,i-1) for i,inc in enumerate(results.iter_visible('increments')): @@ -75,7 +85,7 @@ for filename in options.filenames: shape = [array.shape[0],np.product(array.shape[1:])] vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),deep=True,array_type= vtk.VTK_DOUBLE)) vtk_data[-1].SetName('1_'+x[0].split('/',1)[1]) - rGrid.GetCellData().AddArray(vtk_data[-1]) + grid.GetCellData().AddArray(vtk_data[-1]) else: x = results.get_dataset_location(label) if len(x) == 0: @@ -84,7 +94,7 @@ for filename in options.filenames: shape = [array.shape[0],np.product(array.shape[1:])] vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),deep=True,array_type= vtk.VTK_DOUBLE)) vtk_data[-1].SetName('1_'+x[0].split('/',1)[1]) - rGrid.GetCellData().AddArray(vtk_data[-1]) + grid.GetCellData().AddArray(vtk_data[-1]) results.set_visible('constituents', False) results.set_visible('materialpoints',True) @@ -99,7 +109,7 @@ for filename in options.filenames: shape = [array.shape[0],np.product(array.shape[1:])] vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),deep=True,array_type= vtk.VTK_DOUBLE)) vtk_data[-1].SetName('1_'+x[0].split('/',1)[1]) - rGrid.GetCellData().AddArray(vtk_data[-1]) + grid.GetCellData().AddArray(vtk_data[-1]) else: x = results.get_dataset_location(label) if len(x) == 0: @@ -108,10 +118,10 @@ for filename in options.filenames: shape = [array.shape[0],np.product(array.shape[1:])] vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),deep=True,array_type= vtk.VTK_DOUBLE)) vtk_data[-1].SetName('1_'+x[0].split('/',1)[1]) - rGrid.GetCellData().AddArray(vtk_data[-1]) + grid.GetCellData().AddArray(vtk_data[-1]) - if results.structured: - writer = vtk.vtkXMLRectilinearGridWriter() + writer = vtk.vtkXMLRectilineagridWriter() if results.structured else \ + vtk.vtkXMLUnstructuredGridWriter() results.set_visible('constituents', False) results.set_visible('materialpoints',False) @@ -128,7 +138,6 @@ for filename in options.filenames: writer.SetCompressorTypeToZLib() writer.SetDataModeToBinary() writer.SetFileName(os.path.join(dirname,file_out)) - if results.structured: - writer.SetInputData(rGrid) + writer.SetInputData(grid) writer.Write() diff --git a/processing/post/DADF5_vtk_points.py b/processing/post/DADF5_vtk_points.py new file mode 100755 index 000000000..87c1ad93e --- /dev/null +++ b/processing/post/DADF5_vtk_points.py @@ -0,0 +1,121 @@ +#!/usr/bin/env python3 + +import os +import argparse + +import numpy as np +import vtk +from vtk.util import numpy_support + +import damask + +scriptName = os.path.splitext(os.path.basename(__file__))[0] +scriptID = ' '.join([scriptName,damask.version]) + +# -------------------------------------------------------------------- +# MAIN +# -------------------------------------------------------------------- +parser = argparse.ArgumentParser() + +#ToDo: We need to decide on a way of handling arguments of variable lentght +#https://stackoverflow.com/questions/15459997/passing-integer-lists-to-python + +#parser.add_argument('--version', action='version', version='%(prog)s {}'.format(scriptID)) +parser.add_argument('filenames', nargs='+', + help='DADF5 files') +parser.add_argument('-d','--dir', dest='dir',default='postProc',metavar='string', + help='name of subdirectory relative to the location of the DADF5 file to hold output') +parser.add_argument('--mat', nargs='+', + help='labels for materialpoint',dest='mat') +parser.add_argument('--con', nargs='+', + help='labels for constituent',dest='con') + +options = parser.parse_args() + +if options.mat is None: options.mat=[] +if options.con is None: options.con=[] + +# --- loop over input files ------------------------------------------------------------------------ + +for filename in options.filenames: + results = damask.DADF5(filename) + + Points = vtk.vtkPoints() + Vertices = vtk.vtkCellArray() + for c in results.cell_coordinates(): + pointID = Points.InsertNextPoint(c) + Vertices.InsertNextCell(1) + Vertices.InsertCellPoint(pointID) + + Polydata = vtk.vtkPolyData() + Polydata.SetPoints(Points) + Polydata.SetVerts(Vertices) + Polydata.Modified() + + for i,inc in enumerate(results.iter_visible('increments')): + print('Output step {}/{}'.format(i+1,len(results.increments))) + vtk_data = [] + + results.set_visible('materialpoints',False) + results.set_visible('constituents', True) + for label in options.con: + + for p in results.iter_visible('con_physics'): + if p != 'generic': + for c in results.iter_visible('constituents'): + x = results.get_dataset_location(label) + if len(x) == 0: + continue + array = results.read_dataset(x,0) + shape = [array.shape[0],np.product(array.shape[1:])] + vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),deep=True,array_type= vtk.VTK_DOUBLE)) + vtk_data[-1].SetName('1_'+x[0].split('/',1)[1]) + Polydata.GetCellData().AddArray(vtk_data[-1]) + else: + x = results.get_dataset_location(label) + if len(x) == 0: + continue + array = results.read_dataset(x,0) + shape = [array.shape[0],np.product(array.shape[1:])] + vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),deep=True,array_type= vtk.VTK_DOUBLE)) + vtk_data[-1].SetName('1_'+x[0].split('/',1)[1]) + Polydata.GetCellData().AddArray(vtk_data[-1]) + + results.set_visible('constituents', False) + results.set_visible('materialpoints',True) + for label in options.mat: + for p in results.iter_visible('mat_physics'): + if p != 'generic': + for m in results.iter_visible('materialpoints'): + x = results.get_dataset_location(label) + if len(x) == 0: + continue + array = results.read_dataset(x,0) + shape = [array.shape[0],np.product(array.shape[1:])] + vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),deep=True,array_type= vtk.VTK_DOUBLE)) + vtk_data[-1].SetName('1_'+x[0].split('/',1)[1]) + Polydata.GetCellData().AddArray(vtk_data[-1]) + else: + x = results.get_dataset_location(label) + if len(x) == 0: + continue + array = results.read_dataset(x,0) + shape = [array.shape[0],np.product(array.shape[1:])] + vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),deep=True,array_type= vtk.VTK_DOUBLE)) + vtk_data[-1].SetName('1_'+x[0].split('/',1)[1]) + Polydata.GetCellData().AddArray(vtk_data[-1]) + + writer = vtk.vtkXMLPolyDataWriter() + + + dirname = os.path.abspath(os.path.join(os.path.dirname(filename),options.dir)) + if not os.path.isdir(dirname): + os.mkdir(dirname,0o755) + file_out = '{}_{}.{}'.format(os.path.splitext(os.path.split(filename)[-1])[0],inc,writer.GetDefaultFileExtension()) + + writer.SetCompressorTypeToZLib() + writer.SetDataModeToBinary() + writer.SetFileName(os.path.join(dirname,file_out)) + writer.SetInputData(Polydata) + + writer.Write() diff --git a/processing/post/DADF5toDREAM3D.py b/processing/post/DADF5toDREAM3D.py new file mode 100755 index 000000000..885545297 --- /dev/null +++ b/processing/post/DADF5toDREAM3D.py @@ -0,0 +1,149 @@ +#!/usr/bin/env python3 + +import argparse +import os + +import h5py +import numpy as np + +import damask + +class AttributeManagerNullterm(h5py.AttributeManager): + """ + Attribute management for DREAM.3D hdf5 files. + + String attribute values are stored as fixed-length string with NULLTERM + + References + ---------- + https://stackoverflow.com/questions/38267076 + https://stackoverflow.com/questions/52750232 + + """ + + def create(self, name, data, shape=None, dtype=None): + if isinstance(data,str): + tid = h5py.h5t.C_S1.copy() + tid.set_size(len(data + ' ')) + super().create(name=name,data=data+' ',dtype = h5py.Datatype(tid)) + else: + super().create(name=name,data=data,shape=shape,dtype=dtype) + + +h5py._hl.attrs.AttributeManager = AttributeManagerNullterm # 'Monkey patch' + + +# -------------------------------------------------------------------- +# Crystal structure specifications +# -------------------------------------------------------------------- +Crystal_structures = {'fcc': 1, + 'bcc': 1, + 'hcp': 0, + 'bct': 7, + 'ort': 6} #TODO: is bct Tetragonal low/Tetragonal high? +Phase_types = {'Primary': 0} #further additions to these can be done by looking at 'Create Ensemble Info' filter + + +# -------------------------------------------------------------------- +# MAIN +# -------------------------------------------------------------------- +parser = argparse.ArgumentParser(description='Creating a file for DREAM3D from DAMASK data') +parser.add_argument('filenames',nargs='+',help='HDF5 based output file') +parser.add_argument('--inc',nargs='+',help='Increment for which DREAM3D to be used, eg. 00025',type=int) +parser.add_argument('-d','--dir', dest='dir',default='postProc',metavar='string', + help='name of subdirectory to hold output') + +options = parser.parse_args() + +# -------------------------------------------------------------------- +# loop over input files +for filename in options.filenames: + f = damask.DADF5(filename) #DAMASK output file + count = 0 + for increment in f.increments: + if int(increment[3:]) not in options.inc: + count = count + 1 + continue + + #-------output file creation------------------------------------- + dirname = os.path.abspath(os.path.join(os.path.dirname(filename),options.dir)) + print(dirname) + try: + os.mkdir(dirname) + except FileExistsError: + pass + + o = h5py.File(dirname + '/' + os.path.splitext(filename)[0] + '_{}.dream3D'.format(increment),'w') + #----------------------------------------------------------------- + o.attrs['DADF5toDREAM3D'] = '1.0' + o.attrs['FileVersion'] = '7.0' + #----------------------------------------------------------------- + + + for g in ['DataContainerBundles','Pipeline']: # empty groups (needed) + o.create_group(g) + + data_container_label = 'DataContainers/ImageDataContainer' + cell_data_label = data_container_label + '/CellData' + + + # Phase information of DREAM.3D is constituent ID in DAMASK + o[cell_data_label + '/Phases'] = f.get_constituent_ID().reshape(tuple(f.grid)+(1,)) + # Data quaternions + DAMASK_quaternion = f.read_dataset(f.get_dataset_location('orientation'),0) + DREAM_3D_quaternion = np.empty((np.prod(f.grid),4),dtype=np.float32) + # Convert: DAMASK uses P = -1, DREAM.3D uses P = +1. Also change position of imagninary part + DREAM_3D_quaternion = np.hstack((-DAMASK_quaternion['x'],-DAMASK_quaternion['y'],-DAMASK_quaternion['z'], + DAMASK_quaternion['w'])) + o[cell_data_label + '/Quats'] = DREAM_3D_quaternion.reshape(tuple(f.grid)+(4,)) + + # Attributes to CellData group + o[cell_data_label].attrs['AttributeMatrixType'] = np.array([3],np.uint32) + o[cell_data_label].attrs['TupleDimensions'] = f.grid.astype(np.uint64) + + # Common Attributes for groups in CellData + for group in ['/Phases','/Quats']: + o[cell_data_label + group].attrs['DataArrayVersion'] = np.array([2],np.int32) + o[cell_data_label + group].attrs['Tuple Axis Dimensions'] = 'x={},y={},z={}'.format(*f.grid) + + # phase attributes + o[cell_data_label + '/Phases'].attrs['ComponentDimensions'] = np.array([1],np.uint64) + o[cell_data_label + '/Phases'].attrs['ObjectType'] = 'DataArray' + + # Quats attributes + o[cell_data_label + '/Quats'].attrs['ComponentDimensions'] = np.array([4],np.uint64) + o[cell_data_label + '/Quats'].attrs['ObjectType'] = 'DataArray' + + # Create EnsembleAttributeMatrix + ensemble_label = data_container_label + '/EnsembleAttributeMatrix' + + # Data CrystalStructures + o[ensemble_label + '/CrystalStructures'] = np.uint32(np.array([999,\ + Crystal_structures[f.get_crystal_structure()]])).reshape((2,1)) + o[ensemble_label + '/PhaseTypes'] = np.uint32(np.array([999,Phase_types['Primary']])).reshape((2,1)) # ToDo + + # Attributes Ensemble Matrix + o[ensemble_label].attrs['AttributeMatrixType'] = np.array([11],np.uint32) + o[ensemble_label].attrs['TupleDimensions'] = np.array([2], np.uint64) + + # Attributes for data in Ensemble matrix + for group in ['CrystalStructures','PhaseTypes']: # 'PhaseName' not required MD: But would be nice to take the phase name mapping + o[ensemble_label+'/'+group].attrs['ComponentDimensions'] = np.array([1],np.uint64) + o[ensemble_label+'/'+group].attrs['Tuple Axis Dimensions'] = 'x=2' + o[ensemble_label+'/'+group].attrs['DataArrayVersion'] = np.array([2],np.int32) + o[ensemble_label+'/'+group].attrs['ObjectType'] = 'DataArray' + o[ensemble_label+'/'+group].attrs['TupleDimensions'] = np.array([2],np.uint64) + + + # Create geometry info + geom_label = data_container_label + '/_SIMPL_GEOMETRY' + + o[geom_label + '/DIMENSIONS'] = np.int64(f.grid) + o[geom_label + '/ORIGIN'] = np.float32(np.zeros(3)) + o[geom_label + '/SPACING'] = np.float32(f.size) + + o[geom_label].attrs['GeometryName'] = 'ImageGeometry' + o[geom_label].attrs['GeometryTypeName'] = 'ImageGeometry' + o[geom_label].attrs['GeometryType'] = np.array([0],np.uint32) + o[geom_label].attrs['SpatialDimensionality'] = np.array([3],np.uint32) + o[geom_label].attrs['UnitDimensionality'] = np.array([3],np.uint32) diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index b5e06713f..809c1db98 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -309,6 +309,19 @@ class DADF5(): return path + def get_constituent_ID(self,c=0): + """Pointwise constituent ID.""" + with h5py.File(self.filename,'r') as f: + names = f['/mapping/cellResults/constituent']['Name'][:,c].astype('str') + return np.array([int(n.split('_')[0]) for n in names.tolist()],dtype=np.int32) + + + def get_crystal_structure(self): # ToDo: extension to multi constituents/phase + """Info about the crystal structure.""" + with h5py.File(self.filename,'r') as f: + return f[self.get_dataset_location('orientation')[0]].attrs['Lattice'].astype('str') # np.bytes_ to string + + def read_dataset(self,path,c): """ Dataset for all points/cells. @@ -318,7 +331,7 @@ class DADF5(): with h5py.File(self.filename,'r') as f: shape = (self.Nmaterialpoints,) + np.shape(f[path[0]])[1:] if len(shape) == 1: shape = shape +(1,) - dataset = np.full(shape,np.nan) + dataset = np.full(shape,np.nan,dtype=np.dtype(f[path[0]])) for pa in path: label = pa.split('/')[2] @@ -345,6 +358,20 @@ class DADF5(): return dataset + def cell_coordinates(self): + """Initial coordinates of the cell centers.""" + if self.structured: + delta = self.size/self.grid*0.5 + z, y, x = np.meshgrid(np.linspace(delta[2],self.size[2]-delta[2],self.grid[2]), + np.linspace(delta[1],self.size[1]-delta[1],self.grid[1]), + np.linspace(delta[0],self.size[0]-delta[0],self.grid[0]), + ) + return np.concatenate((x[:,:,:,None],y[:,:,:,None],y[:,:,:,None]),axis = 3).reshape([np.product(self.grid),3]) + else: + with h5py.File(self.filename,'r') as f: + return f['geometry/x_c'][()] + + def add_Cauchy(self,P='P',F='F'): """ Adds Cauchy stress calculated from 1st Piola-Kirchhoff stress and deformation gradient. diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 6649de542..32098b109 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -381,6 +381,7 @@ subroutine CPFEM_results(inc,time) call results_addIncrement(inc,time) call constitutive_results call crystallite_results + call homogenization_results call results_removeLink('current') ! ToDo: put this into closeJobFile call results_closeJobFile #endif diff --git a/src/DAMASK_abaqus.f b/src/DAMASK_abaqus.f index 514307fe8..7a43f688e 100644 --- a/src/DAMASK_abaqus.f +++ b/src/DAMASK_abaqus.f @@ -313,11 +313,10 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,& call CPFEM_general(computationMode,usePingPong,dfgrd0,dfgrd1,temperature,dtime,noel,npt,stress_h,ddsdde_h) -! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13 -! straight: 11, 22, 33, 12, 23, 13 -! ABAQUS explicit: 11, 22, 33, 12, 23, 13 -! ABAQUS implicit: 11, 22, 33, 12, 13, 23 -! ABAQUS implicit: 11, 22, 33, 12 +! DAMASK: 11, 22, 33, 12, 23, 13 +! ABAQUS explicit: 11, 22, 33, 12, 23, 13 +! ABAQUS implicit: 11, 22, 33, 12, 13, 23 +! ABAQUS implicit: 11, 22, 33, 12 ddsdde = ddsdde_h(1:ntens,1:ntens) stress = stress_h(1:ntens) diff --git a/src/DAMASK_interface.f90 b/src/DAMASK_interface.f90 index 7febcef13..ee8f220ae 100644 --- a/src/DAMASK_interface.f90 +++ b/src/DAMASK_interface.f90 @@ -13,7 +13,7 @@ #define INTEL_MIN 1600 #define PETSC_MAJOR 3 #define PETSC_MINOR_MIN 10 -#define PETSC_MINOR_MAX 11 +#define PETSC_MINOR_MAX 12 module DAMASK_interface use, intrinsic :: iso_fortran_env diff --git a/src/DAMASK_marc.f90 b/src/DAMASK_marc.f90 index 3f527daf5..09a1c83c8 100644 --- a/src/DAMASK_marc.f90 +++ b/src/DAMASK_marc.f90 @@ -265,8 +265,8 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & call debug_reset() ! resets debugging outdatedFFN1 = .false. cycleCounter = cycleCounter + 1 - mesh_cellnode = mesh_build_cellnodes() ! update cell node coordinates - call mesh_build_ipCoordinates() ! update ip coordinates + !mesh_cellnode = mesh_build_cellnodes() ! update cell node coordinates + !call mesh_build_ipCoordinates() ! update ip coordinates endif if (outdatedByNewInc) then computationMode = ior(computationMode,CPFEM_AGERESULTS) ! calc and age results @@ -315,9 +315,6 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & lastLovl = lovl ! record lovl call CPFEM_general(computationMode,usePingPong,ffn,ffn1,t(1),timinc,m(1),nn,stress,ddsdde) -! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13 -! Marc: 11, 22, 33, 12, 23, 13 -! Marc: 11, 22, 33, 12 d = ddsdde(1:ngens,1:ngens) s = stress(1:ndi+nshear) diff --git a/src/HDF5_utilities.f90 b/src/HDF5_utilities.f90 index 68281a6a0..e4819431e 100644 --- a/src/HDF5_utilities.f90 +++ b/src/HDF5_utilities.f90 @@ -112,14 +112,10 @@ subroutine HDF5_utilities_init call h5open_f(hdferr) if (hdferr < 0) call IO_error(1,ext_msg='HDF5_Utilities_init: h5open_f') -#ifndef Marc4DAMASK -! This test should ensure that integer size matches. For some reasons, the HDF5 libraries -! that come with MSC.Marc>=2019 seem to be of 4byte even though it is a 8byte Marc version call h5tget_size_f(H5T_NATIVE_INTEGER,typeSize, hdferr) if (hdferr < 0) call IO_error(1,ext_msg='HDF5_Utilities_init: h5tget_size_f (int)') if (int(bit_size(0),SIZE_T)/=typeSize*8) & call IO_error(0,ext_msg='Default integer size does not match H5T_NATIVE_INTEGER') -#endif call h5tget_size_f(H5T_NATIVE_DOUBLE,typeSize, hdferr) if (hdferr < 0) call IO_error(1,ext_msg='HDF5_Utilities_init: h5tget_size_f (double)') diff --git a/src/IO.f90 b/src/IO.f90 index 066ebbe8c..5d96e67d2 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -44,7 +44,6 @@ module IO IO_extractValue, & IO_countDataLines #elif defined(Marc4DAMASK) - IO_skipChunks, & IO_fixedNoEFloatValue, & IO_fixedIntValue, & IO_countNumericalDataLines @@ -189,17 +188,17 @@ integer function IO_open_binary(fileName,mode) m = 'r' endif - if (m == 'w') then - open(newunit=IO_open_binary, file=trim(fileName),& - status='replace',access='stream',action='write',iostat=ierr) - if (ierr /= 0) call IO_error(100,ext_msg='could not open file (w): '//trim(fileName)) - elseif(m == 'r') then - open(newunit=IO_open_binary, file=trim(fileName),& - status='old', access='stream',action='read', iostat=ierr) - if (ierr /= 0) call IO_error(100,ext_msg='could not open file (r): '//trim(fileName)) - else - call IO_error(100,ext_msg='unknown access mode: '//m) - endif + if (m == 'w') then + open(newunit=IO_open_binary, file=trim(fileName),& + status='replace',access='stream',action='write',iostat=ierr) + if (ierr /= 0) call IO_error(100,ext_msg='could not open file (w): '//trim(fileName)) + elseif(m == 'r') then + open(newunit=IO_open_binary, file=trim(fileName),& + status='old', access='stream',action='read', iostat=ierr) + if (ierr /= 0) call IO_error(100,ext_msg='could not open file (r): '//trim(fileName)) + else + call IO_error(100,ext_msg='unknown access mode: '//m) + endif end function IO_open_binary @@ -403,7 +402,7 @@ pure function IO_stringPos(string) left = right + verify(string(right+1:),SEP) right = left + scan(string(left:),SEP) - 2 if ( string(left:left) == '#' ) exit - IO_stringPos = [IO_stringPos,left, right] + IO_stringPos = [IO_stringPos,left,right] IO_stringPos(1) = IO_stringPos(1)+1 endOfString: if (right < left) then IO_stringPos(IO_stringPos(1)*2+1) = len_trim(string) @@ -1023,27 +1022,6 @@ integer function IO_countNumericalDataLines(fileUnit) backspace(fileUnit) end function IO_countNumericalDataLines - - -!-------------------------------------------------------------------------------------------------- -!> @brief reads file to skip (at least) N chunks (may be over multiple lines) -!-------------------------------------------------------------------------------------------------- -subroutine IO_skipChunks(fileUnit,N) - - integer, intent(in) :: fileUnit, & !< file handle - N !< minimum number of chunks to skip - - integer :: remainingChunks - character(len=65536) :: line - - line = '' - remainingChunks = N - - do while (trim(line) /= IO_EOF .and. remainingChunks > 0) - line = IO_read(fileUnit) - remainingChunks = remainingChunks - (size(IO_stringPos(line))-1)/2 - enddo -end subroutine IO_skipChunks #endif @@ -1072,8 +1050,8 @@ integer function IO_countContinuousIntValues(fileUnit) if (chunkPos(1) < 1) then ! empty line exit elseif (IO_lc(IO_stringValue(line,chunkPos,2)) == 'to' ) then ! found range indicator - IO_countContinuousIntValues = 1 + abs( IO_intValue(line,chunkPos,3) & - - IO_intValue(line,chunkPos,1)) + IO_countContinuousIntValues = 1 + abs( IO_intValue(line,chunkPos,3) & + -IO_intValue(line,chunkPos,1)) exit ! only one single range indicator allowed else IO_countContinuousIntValues = IO_countContinuousIntValues+chunkPos(1)-1 ! add line's count when assuming 'c' diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index ac4aa009a..79a385818 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -14,11 +14,11 @@ #include "Lambert.f90" #include "rotations.f90" #include "FEsolving.f90" -#include "geometry_plastic_nonlocal.f90" #include "element.f90" #include "mesh_base.f90" #include "HDF5_utilities.f90" #include "results.f90" +#include "geometry_plastic_nonlocal.f90" #include "discretization.f90" #ifdef Abaqus #include "mesh_abaqus.f90" diff --git a/src/crystallite.f90 b/src/crystallite.f90 index a31b22db9..f3459c0d7 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -2090,7 +2090,7 @@ end subroutine setConvergenceFlag !-------------------------------------------------------------------------------------------------- -!> @brief Standard forwarding of state as state = state0 + dotState * (delta t) +!> @brief Standard forwarding of state as state = state0 + dotState * (delta t) comment seems wrong! !-------------------------------------------------------------------------------------------------- subroutine update_stress(timeFraction) diff --git a/src/element.f90 b/src/element.f90 index d62b4fd93..02f5fb762 100644 --- a/src/element.f90 +++ b/src/element.f90 @@ -146,11 +146,11 @@ module element 8 & ! 3D 8node ] !< number of cell nodes in a specific cell type - ! *** FE_ipNeighbor *** - ! is a list of the neighborhood of each IP. + ! *** IPneighbor *** + ! list of the neighborhood of each IP. ! It is sorted in (local) +x,-x, +y,-y, +z,-z direction. - ! Positive integers denote an intra-FE IP identifier. - ! Negative integers denote the interface behind which the neighboring (extra-FE) IP will be located. + ! Positive integers denote an intra-element IP identifier. + ! Negative integers denote the interface behind which the neighboring (extra-element) IP will be located. integer, dimension(nIPneighbor(cellType(1)),nIP(1)), parameter, private :: IPneighbor1 = & reshape([& @@ -256,10 +256,6 @@ module element -3,26,-4,24,-6,18 & ],[nIPneighbor(cellType(10)),nIP(10)]) -! MD: probably not needed END -! -------------------------------------------------------------------------------------------------- - - integer, dimension(nNode(1),NcellNode(geomType(1))), parameter :: cellNodeParentNodeWeights1 = & reshape([& @@ -757,7 +753,7 @@ subroutine tElement_init(self,elemType) self%cell = CELL10 end select - self%NcellNodesPerCell = NCELLNODEPERCELL(self%cellType) + self%NcellnodesPerCell = NCELLNODEPERCELL(self%cellType) select case(self%cellType) case(1) diff --git a/src/geometry_plastic_nonlocal.f90 b/src/geometry_plastic_nonlocal.f90 index 37909a4a5..88634c245 100644 --- a/src/geometry_plastic_nonlocal.f90 +++ b/src/geometry_plastic_nonlocal.f90 @@ -7,6 +7,7 @@ !-------------------------------------------------------------------------------------------------- module geometry_plastic_nonlocal use prec + use results implicit none private @@ -32,6 +33,7 @@ module geometry_plastic_nonlocal geometry_plastic_nonlocal_setIPvolume, & geometry_plastic_nonlocal_setIParea, & geometry_plastic_nonlocal_setIPareaNormal, & + geometry_plastic_nonlocal_results, & geometry_plastic_nonlocal_disable contains @@ -112,4 +114,45 @@ subroutine geometry_plastic_nonlocal_disable end subroutine geometry_plastic_nonlocal_disable + +!--------------------------------------------------------------------------------------------------- +!> @brief Writes geometry data to results file +!--------------------------------------------------------------------------------------------------- +subroutine geometry_plastic_nonlocal_results + + integer, dimension(:), allocatable :: shp + +#if defined(DAMASK_HDF5) + call results_openJobFile + + writeVolume: block + real(pReal), dimension(:), allocatable :: temp + shp = shape(geometry_plastic_nonlocal_IPvolume0) + temp = reshape(geometry_plastic_nonlocal_IPvolume0,[shp(1)*shp(2)]) + call results_writeDataset('geometry',temp,'v_0',& + 'initial cell volume','m³') + end block writeVolume + + writeAreas: block + real(pReal), dimension(:,:), allocatable :: temp + shp = shape(geometry_plastic_nonlocal_IParea0) + temp = reshape(geometry_plastic_nonlocal_IParea0,[shp(1),shp(2)*shp(3)]) + call results_writeDataset('geometry',temp,'a_0',& + 'initial cell face area','m²') + end block writeAreas + + writeNormals: block + real(pReal), dimension(:,:,:), allocatable :: temp + shp = shape(geometry_plastic_nonlocal_IPareaNormal0) + temp = reshape(geometry_plastic_nonlocal_IPareaNormal0,[shp(1),shp(2),shp(3)*shp(4)]) + call results_writeDataset('geometry',temp,'n_0',& + 'initial cell face normals','-',transposed=.false.) + end block writeNormals + + + call results_closeJobFile +#endif + +end subroutine geometry_plastic_nonlocal_results + end module geometry_plastic_nonlocal diff --git a/src/lattice.f90 b/src/lattice.f90 index c6d559790..a09355de3 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -228,92 +228,92 @@ module lattice real(pReal), dimension(4+4,LATTICE_HEX_NSLIP), parameter :: & LATTICE_HEX_SYSTEMSLIP = reshape(real([& ! Slip direction Plane normal - ! Basal systems <11.0>{00.1} (independent of c/a-ratio, Bravais notation (4 coordinate base)) + ! Basal systems <-1-1.0>{00.1} (independent of c/a-ratio, Bravais notation (4 coordinate base)) 2, -1, -1, 0, 0, 0, 0, 1, & -1, 2, -1, 0, 0, 0, 0, 1, & -1, -1, 2, 0, 0, 0, 0, 1, & - ! 1st type prismatic systems <11.0>{10.0} (independent of c/a-ratio) + ! 1st type prismatic systems <-1-1.0>{1-1.0} (independent of c/a-ratio) 2, -1, -1, 0, 0, 1, -1, 0, & -1, 2, -1, 0, -1, 0, 1, 0, & -1, -1, 2, 0, 1, -1, 0, 0, & - ! 2nd type prismatic systems <10.0>{11.0} -- a slip; plane normals independent of c/a-ratio - 0, 1, -1, 0, 2, -1, -1, 0, & - -1, 0, 1, 0, -1, 2, -1, 0, & - 1, -1, 0, 0, -1, -1, 2, 0, & - ! 1st type 1st order pyramidal systems <11.0>{-11.1} -- plane normals depend on the c/a-ratio - 2, -1, -1, 0, 0, 1, -1, 1, & - -1, 2, -1, 0, -1, 0, 1, 1, & - -1, -1, 2, 0, 1, -1, 0, 1, & - 1, 1, -2, 0, -1, 1, 0, 1, & - -2, 1, 1, 0, 0, -1, 1, 1, & - 1, -2, 1, 0, 1, 0, -1, 1, & + ! 2nd type prismatic systems <-11.0>{11.0} -- a slip; plane normals independent of c/a-ratio + -1, 1, 0, 0, 1, 1, -2, 0, & + 0, -1, 1, 0, -2, 1, 1, 0, & + 1, 0, -1, 0, 1, -2, 1, 0, & + ! 1st type 1st order pyramidal systems <-1-1.0>{-11.1} -- plane normals depend on the c/a-ratio + -1, 2, -1, 0, 1, 0, -1, 1, & + -2, 1, 1, 0, 0, 1, -1, 1, & + -1, -1, 2, 0, -1, 1, 0, 1, & + 1, -2, 1, 0, -1, 0, 1, 1, & + 2, -1, -1, 0, 0, -1, 1, 1, & + 1, 1, -2, 0, 1, -1, 0, 1, & ! pyramidal system: c+a slip <11.3>{-10.1} -- plane normals depend on the c/a-ratio - 2, -1, -1, 3, -1, 1, 0, 1, & - 1, -2, 1, 3, -1, 1, 0, 1, & - -1, -1, 2, 3, 1, 0, -1, 1, & -2, 1, 1, 3, 1, 0, -1, 1, & - -1, 2, -1, 3, 0, -1, 1, 1, & - 1, 1, -2, 3, 0, -1, 1, 1, & - -2, 1, 1, 3, 1, -1, 0, 1, & - -1, 2, -1, 3, 1, -1, 0, 1, & - 1, 1, -2, 3, -1, 0, 1, 1, & - 2, -1, -1, 3, -1, 0, 1, 1, & - 1, -2, 1, 3, 0, 1, -1, 1, & + -1, -1, 2, 3, 1, 0, -1, 1, & -1, -1, 2, 3, 0, 1, -1, 1, & + 1, -2, 1, 3, 0, 1, -1, 1, & + 1, -2, 1, 3, -1, 1, 0, 1, & + 2, -1, -1, 3, -1, 1, 0, 1, & + 2, -1, -1, 3, -1, 0, 1, 1, & + 1, 1, -2, 3, -1, 0, 1, 1, & + 1, 1, -2, 3, 0, -1, 1, 1, & + -1, 2, -1, 3, 0, -1, 1, 1, & + -1, 2, -1, 3, 1, -1, 0, 1, & + -2, 1, 1, 3, 1, -1, 0, 1, & ! pyramidal system: c+a slip <11.3>{-1-1.2} -- as for hexagonal ice (Castelnau et al. 1996, similar to twin system found below) - 2, -1, -1, 3, -2, 1, 1, 2, & ! sorted according to similar twin system - -1, 2, -1, 3, 1, -2, 1, 2, & ! <11.3>{-1-1.2} shear = 2((c/a)^2-2)/(3 c/a) - -1, -1, 2, 3, 1, 1, -2, 2, & - -2, 1, 1, 3, 2, -1, -1, 2, & + -1, -1, 2, 3, 1, 1, -2, 2, & ! <11.3>{-1-1.2} shear = 2((c/a)^2-2)/(3 c/a) 1, -2, 1, 3, -1, 2, -1, 2, & - 1, 1, -2, 3, -1, -1, 2, 2 & - ],pReal),shape(LATTICE_HEX_SYSTEMSLIP)) !< slip systems for hex sorted by A. Alankar & P. Eisenlohr + 2, -1, -1, 3, -2, 1, 1, 2, & + 1, 1, -2, 3, -1, -1, 2, 2, & + -1, 2, -1, 3, 1, -2, 1, 2, & + -2, 1, 1, 3, 2, -1, -1, 2 & + ],pReal),shape(LATTICE_HEX_SYSTEMSLIP)) !< slip systems for hex, sorted by P. Eisenlohr CCW around starting next to a_1 axis character(len=*), dimension(6), parameter :: LATTICE_HEX_SLIPFAMILY_NAME = & - ['<1 1 . 1>{0 0 . 1} ', & - '<1 1 . 1>{1 0 . 0} ', & - '<1 0 . 0>{1 1 . 0} ', & - '<1 1 . 0>{-1 1 . 1} ', & - '<1 1 . 3>{-1 0 . 1} ', & - '<1 1 . 3>{-1 -1 . 2}'] + ['< 1 1 . 0>{ 0 0 . 1}', & + '< 1 1 . 0>{ 1 0 . 0}', & + '<-1 1 . 0>{ 1 1 . 0}', & + '< 1 1 . 0>{ 1 -1 . 1}', & + '< 1 1 . 3>{-1 0 . 1}', & + '< 1 1 . 3>{-1 -1 . 2}'] real(pReal), dimension(4+4,LATTICE_HEX_NTWIN), parameter :: & LATTICE_HEX_SYSTEMTWIN = reshape(real([& - ! Compression or Tension =f(twinning shear=f(c/a)) for each metal ! (according to Yoo 1981) - 1, -1, 0, 1, -1, 1, 0, 2, & ! <-10.1>{10.2} shear = (3-(c/a)^2)/(sqrt(3) c/a) - -1, 0, 1, 1, 1, 0, -1, 2, & + ! Compression or Tension = f(twinning shear=f(c/a)) for each metal ! (according to Yoo 1981) + -1, 0, 1, 1, 1, 0, -1, 2, & ! <-10.1>{10.2} shear = (3-(c/a)^2)/(sqrt(3) c/a) + 0, -1, 1, 1, 0, 1, -1, 2, & + 1, -1, 0, 1, -1, 1, 0, 2, & + 1, 0, -1, 1, -1, 0, 1, 2, & 0, 1, -1, 1, 0, -1, 1, 2, & -1, 1, 0, 1, 1, -1, 0, 2, & - 1, 0, -1, 1, -1, 0, 1, 2, & - 0, -1, 1, 1, 0, 1, -1, 2, & ! - 2, -1, -1, 6, -2, 1, 1, 1, & ! <11.6>{-1-1.1} shear = 1/(c/a) - -1, 2, -1, 6, 1, -2, 1, 1, & - -1, -1, 2, 6, 1, 1, -2, 1, & - -2, 1, 1, 6, 2, -1, -1, 1, & + -1, -1, 2, 6, 1, 1, -2, 1, & ! <11.6>{-1-1.1} shear = 1/(c/a) 1, -2, 1, 6, -1, 2, -1, 1, & + 2, -1, -1, 6, -2, 1, 1, 1, & 1, 1, -2, 6, -1, -1, 2, 1, & + -1, 2, -1, 6, 1, -2, 1, 1, & + -2, 1, 1, 6, 2, -1, -1, 1, & ! - -1, 1, 0, -2, -1, 1, 0, 1, & !! <10.-2>{10.1} shear = (4(c/a)^2-9)/(4 sqrt(3) c/a) - 1, 0, -1, -2, 1, 0, -1, 1, & + 1, 0, -1, -2, 1, 0, -1, 1, & ! <10.-2>{10.1} shear = (4(c/a)^2-9)/(4 sqrt(3) c/a) + 0, 1, -1, -2, 0, 1, -1, 1, & + -1, 1, 0, -2, -1, 1, 0, 1, & + -1, 0, 1, -2, -1, 0, 1, 1, & 0, -1, 1, -2, 0, -1, 1, 1, & 1, -1, 0, -2, 1, -1, 0, 1, & - -1, 0, 1, -2, -1, 0, 1, 1, & - 0, 1, -1, -2, 0, 1, -1, 1, & ! - 2, -1, -1, -3, 2, -1, -1, 2, & ! <11.-3>{11.2} shear = 2((c/a)^2-2)/(3 c/a) + 1, 1, -2, -3, 1, 1, -2, 2, & ! <11.-3>{11.2} shear = 2((c/a)^2-2)/(3 c/a) -1, 2, -1, -3, -1, 2, -1, 2, & - -1, -1, 2, -3, -1, -1, 2, 2, & -2, 1, 1, -3, -2, 1, 1, 2, & + -1, -1, 2, -3, -1, -1, 2, 2, & 1, -2, 1, -3, 1, -2, 1, 2, & - 1, 1, -2, -3, 1, 1, -2, 2 & - ],pReal),shape(LATTICE_HEX_SYSTEMTWIN)) !< twin systems for hex, order follows Prof. Tom Bieler's scheme + 2, -1, -1, -3, 2, -1, -1, 2 & + ],pReal),shape(LATTICE_HEX_SYSTEMTWIN)) !< twin systems for hex, sorted by P. Eisenlohr CCW around starting next to a_1 axis character(len=*), dimension(4), parameter :: LATTICE_HEX_TWINFAMILY_NAME = & - ['<-1 0 . 1>{1 0 . 2} ', & - '<1 1 . 6>{-1 -1 . 1}', & - '<1 0 . -2>{1 0 . 1} ', & - '<1 1 . -3>{1 1 . 2} '] + ['<-1 0 . 1>{ 1 0 . 2}', & + '< 1 1 . 6>{-1 -1 . 1}', & + '< 1 0 . -2>{ 1 0 . 1}', & + '< 1 1 . -3>{ 1 1 . 2}'] real(pReal), dimension(4+4,LATTICE_HEX_NCLEAVAGE), parameter :: & LATTICE_HEX_SYSTEMCLEAVAGE = reshape(real([& diff --git a/src/math.f90 b/src/math.f90 index 26e63388c..5a3a5f467 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -1307,10 +1307,10 @@ real(pReal) pure function math_volTetrahedron(v1,v2,v3,v4) real(pReal), dimension (3,3) :: m m(1:3,1) = v1-v2 - m(1:3,2) = v2-v3 - m(1:3,3) = v3-v4 + m(1:3,2) = v1-v3 + m(1:3,3) = v1-v4 - math_volTetrahedron = math_det33(m)/6.0_pReal + math_volTetrahedron = abs(math_det33(m))/6.0_pReal end function math_volTetrahedron @@ -1404,6 +1404,7 @@ subroutine unitTest integer, dimension(5) :: range_out_ = [1,2,3,4,5] real(pReal) :: det + real(pReal), dimension(3) :: v3_1,v3_2,v3_3,v3_4 real(pReal), dimension(6) :: v6 real(pReal), dimension(9) :: v9 real(pReal), dimension(3,3) :: t33,t33_2 @@ -1452,6 +1453,15 @@ subroutine unitTest if(any(dNeq0(math_6toSym33(v6) - math_symmetric33(math_6toSym33(v6))))) & call IO_error(401,ext_msg='math_symmetric33') + call random_number(v3_1) + call random_number(v3_2) + call random_number(v3_3) + call random_number(v3_4) + + if(dNeq(abs(dot_product(math_cross(v3_1-v3_4,v3_2-v3_4),v3_3-v3_4))/6.0, & + math_volTetrahedron(v3_1,v3_2,v3_3,v3_4),tol=1.0e-12_pReal)) & + call IO_error(401,ext_msg='math_volTetrahedron') + call random_number(t33) if(dNeq(math_det33(math_symmetric33(t33)),math_detSym33(math_symmetric33(t33)),tol=1.0e-12_pReal)) & call IO_error(401,ext_msg='math_det33/math_detSym33') diff --git a/src/mesh_FEM.f90 b/src/mesh_FEM.f90 index 2fd63c0b3..4b4ff18aa 100644 --- a/src/mesh_FEM.f90 +++ b/src/mesh_FEM.f90 @@ -28,8 +28,7 @@ module mesh integer, public, protected :: & mesh_Nboundaries, & mesh_NcpElems, & !< total number of CP elements in mesh - mesh_NcpElemsGlobal, & - mesh_Nnodes !< total number of nodes in mesh + mesh_NcpElemsGlobal !!!! BEGIN DEPRECATED !!!!! integer, public, protected :: & @@ -69,8 +68,7 @@ module mesh public :: & mesh_init, & mesh_FEM_build_ipVolumes, & - mesh_FEM_build_ipCoordinates, & - mesh_cellCenterCoordinates + mesh_FEM_build_ipCoordinates contains @@ -107,7 +105,8 @@ subroutine mesh_init integer, parameter :: FILEUNIT = 222 integer :: j integer, allocatable, dimension(:) :: chunkPos - integer :: dimPlex + integer :: dimPlex, & + mesh_Nnodes !< total number of nodes in mesh integer, parameter :: & mesh_ElemType=1 !< Element type of the mesh (only support homogeneous meshes) character(len=512) :: & @@ -221,9 +220,6 @@ subroutine mesh_init call theMesh%init(dimplex,integrationOrder,mesh_node0) call theMesh%setNelems(mesh_NcpElems) - theMesh%homogenizationAt = mesh_element(3,:) - theMesh%microstructureAt = mesh_element(4,:) - call discretization_init(mesh_element(3,:),mesh_element(4,:),& reshape(mesh_ipCoordinates,[3,mesh_maxNips*mesh_NcpElems]), & mesh_node0) @@ -231,26 +227,8 @@ subroutine mesh_init end subroutine mesh_init -!-------------------------------------------------------------------------------------------------- -!> @brief Calculates cell center coordinates. -!-------------------------------------------------------------------------------------------------- -pure function mesh_cellCenterCoordinates(ip,el) - - integer, intent(in) :: el, & !< element number - ip !< integration point number - real(pReal), dimension(3) :: mesh_cellCenterCoordinates !< x,y,z coordinates of the cell center of the requested IP cell - -end function mesh_cellCenterCoordinates - - !-------------------------------------------------------------------------------------------------- !> @brief Calculates IP volume. Allocates global array 'mesh_ipVolume' -!> @details The IP volume is calculated differently depending on the cell type. -!> 2D cells assume an element depth of one in order to calculate the volume. -!> For the hexahedral cell we subdivide the cell into subvolumes of pyramidal -!> shape with a cell face as basis and the central ip at the tip. This subvolume is -!> calculated as an average of four tetrahedals with three corners on the cell face -!> and one corner at the central ip. !-------------------------------------------------------------------------------------------------- subroutine mesh_FEM_build_ipVolumes(dimPlex) diff --git a/src/mesh_base.f90 b/src/mesh_base.f90 index 102f53e9e..dab7059ee 100644 --- a/src/mesh_base.f90 +++ b/src/mesh_base.f90 @@ -1,4 +1,3 @@ - !-------------------------------------------------------------------------------------------------- !> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH @@ -8,14 +7,13 @@ !-------------------------------------------------------------------------------------------------- module mesh_base - use, intrinsic :: iso_c_binding use prec use element implicit none !--------------------------------------------------------------------------------------------------- -!> Properties of a the whole mesh (consisting of one type of elements) +!> Properties of a whole mesh (consisting of one type of elements) !--------------------------------------------------------------------------------------------------- type, public :: tMesh type(tElement) :: & @@ -33,11 +31,7 @@ module mesh_base elemType, & Ncells, & nIPneighbors, & - NcellNodes, & - maxElemsPerNode - integer(pInt), dimension(:), allocatable, public :: & - homogenizationAt, & - microstructureAt + NcellNodes integer(pInt), dimension(:,:), allocatable, public :: & connectivity contains @@ -47,6 +41,7 @@ module mesh_base end type tMesh contains + subroutine tMesh_base_init(self,meshType,elemType,nodes) class(tMesh) :: self diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 4a5f17674..66906207a 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -23,118 +23,29 @@ module mesh implicit none private - + + type tCellNodeDefinition + integer, dimension(:,:), allocatable :: parents + integer, dimension(:,:), allocatable :: weights + end type tCellNodeDefinition + + type(tCellNodeDefinition), dimension(:), allocatable :: cellNodeDefinition + real(pReal), public, protected :: & - mesh_unitlength !< physical length of one unit in mesh + mesh_unitlength !< physical length of one unit in mesh + + integer, dimension(:,:), allocatable, target :: & + mesh_mapFEtoCPelem, & !< [sorted FEid, corresponding CPid] + mesh_mapFEtoCPnode !< [sorted FEid, corresponding CPid] !-------------------------------------------------------------------------------------------------- -! public variables (DEPRECATED) - - real(pReal), dimension(:,:,:), allocatable, public :: & - mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!) - - real(pReal), dimension(:,:), allocatable, public :: & - mesh_cellnode !< cell node x,y,z coordinates (after deformation! ONLY FOR MARC!!!) +! DEPRECATED + real(pReal), dimension(:,:,:), allocatable, public :: & + mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!) !-------------------------------------------------------------------------------------------------- - integer, dimension(:,:), allocatable :: & - mesh_element - - integer, dimension(:,:,:,:), allocatable :: & - mesh_ipNeighborhood !< 6 or less neighboring IPs as [element_num, IP_index, neighbor_index that points to me] - - real(pReal), dimension(:,:), allocatable :: & - mesh_node, & !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!! - mesh_ipVolume, & !< volume associated with IP (initially!) - mesh_node0 !< node x,y,z coordinates (initially!) - - real(pReal), dimension(:,:,:), allocatable:: & - mesh_ipArea !< area of interface to neighboring IP (initially!) - real(pReal),dimension(:,:,:,:), allocatable :: & - mesh_ipAreaNormal !< area normal of interface to neighboring IP (initially!) - -! -------------------------------------------------------------------------------------------------- - - type(tMesh) :: theMesh - - - integer:: & - mesh_Ncellnodes, & !< total number of cell nodes in mesh (including duplicates) - mesh_elemType, & !< Element type of the mesh (only support homogeneous meshes) - mesh_Nnodes, & !< total number of nodes in mesh - mesh_Ncells, & !< total number of cells in mesh - mesh_maxNsharedElems !< max number of CP elements sharing a node - - -integer, dimension(:,:), allocatable :: & - mesh_cellnodeParent !< cellnode's parent element ID, cellnode's intra-element ID - - integer,dimension(:,:,:), allocatable :: & - mesh_cell2, & !< cell connectivity for each element,ip/cell - mesh_cell !< cell connectivity for each element,ip/cell - -! These definitions should actually reside in the FE-solver specific part (different for MARC/ABAQUS) -! Hence, I suggest to prefix with "FE_" - - integer, parameter :: & - FE_Ngeomtypes = 10, & - FE_Ncelltypes = 4, & - FE_maxNcellnodesPerCell = 8, & - FE_maxNcellnodesPerCellface = 4 - - integer, dimension(FE_Ngeomtypes), parameter :: FE_NmatchingNodes = & !< number of nodes that are needed for face matching in a specific type of element geometry - int([ & - 3, & ! element 6 (2D 3node 1ip) - 3, & ! element 125 (2D 6node 3ip) - 4, & ! element 11 (2D 4node 4ip) - 4, & ! element 27 (2D 8node 9ip) - 4, & ! element 134 (3D 4node 1ip) - 4, & ! element 127 (3D 10node 4ip) - 6, & ! element 136 (3D 6node 6ip) - 8, & ! element 117 (3D 8node 1ip) - 8, & ! element 7 (3D 8node 8ip) - 8 & ! element 21 (3D 20node 27ip) - ],pInt) - - - integer, dimension(FE_Ncelltypes), parameter :: FE_NcellnodesPerCellface = & !< number of cell nodes per cell face in a specific cell type - int([& - 2, & ! (2D 3node) - 2, & ! (2D 4node) - 3, & ! (3D 4node) - 4 & ! (3D 8node) - ],pInt) - - integer, dimension(FE_Ncelltypes), parameter :: FE_NipNeighbors = & !< number of ip neighbors / cell faces in a specific cell type - int([& - 3, & ! (2D 3node) - 4, & ! (2D 4node) - 4, & ! (3D 4node) - 6 & ! (3D 8node) - ],pInt) - - - integer :: & - mesh_NelemSets - character(len=64), dimension(:), allocatable :: & - mesh_nameElemSet - - integer, dimension(:,:), allocatable :: & - mesh_mapElemSet !< list of elements in elementSet - integer, dimension(:,:), allocatable, target :: & - mesh_mapFEtoCPelem, & !< [sorted FEid, corresponding CPid] - mesh_mapFEtoCPnode !< [sorted FEid, corresponding CPid] - - integer, dimension(:,:,:,:), allocatable :: & - mesh_ipNeighborhood2 !< 6 or less neighboring IPs as [element_num, IP_index, neighbor_index that points to me] - - integer, dimension(:), allocatable :: & - Marc_matNumber !< array of material numbers for hypoelastic material (Marc only) - public :: & mesh_init, & - mesh_build_cellnodes, & - mesh_build_ipCoordinates, & mesh_FEasCP @@ -148,241 +59,305 @@ subroutine mesh_init(ip,el) integer, intent(in) :: el, ip - integer, parameter :: FILEUNIT = 222 - integer :: j, fileFormatVersion, elemType, & - mesh_maxNelemInSet, & - mesh_nElems, & - hypoelasticTableStyle, & - initialcondTableStyle - logical :: myDebug - + real(pReal), dimension(:,:), allocatable :: & + node0_elem, & !< node x,y,z coordinates (initially!) + node0_cell + type(tElement) :: elem + + integer :: nElems + integer, dimension(:), allocatable :: & + microstructureAt, & + homogenizationAt + integer:: & + Nnodes !< total number of nodes in mesh + + real(pReal), dimension(:,:), allocatable :: & + ip_reshaped + integer,dimension(:,:,:), allocatable :: & + connectivity_cell !< cell connectivity for each element,ip/cell + integer, dimension(:,:), allocatable :: & + connectivity_elem + real(pReal), dimension(:,:,:,:),allocatable :: & + unscaledNormals + write(6,'(/,a)') ' <<<+- mesh init -+>>>' - mesh_unitlength = numerics_unitlength ! set physical extent of a length unit in mesh + mesh_unitlength = numerics_unitlength ! set physical extent of a length unit in mesh + + call inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogenizationAt) + nElems = size(connectivity_elem,2) + + if (debug_e < 1 .or. debug_e > nElems) call IO_error(602,ext_msg='element') + if (debug_i < 1 .or. debug_i > elem%nIPs) call IO_error(602,ext_msg='IP') - myDebug = (iand(debug_level(debug_mesh),debug_levelBasic) /= 0) + FEsolving_execElem = [ 1,nElems ] ! parallel loop bounds set to comprise all DAMASK elements + allocate(FEsolving_execIP(2,nElems), source=1) ! parallel loop bounds set to comprise from first IP... + FEsolving_execIP(2,:) = elem%nIPs - call IO_open_inputFile(FILEUNIT,modelName) - fileFormatVersion = mesh_marc_get_fileFormat(FILEUNIT) - if (myDebug) write(6,'(a)') ' Got input file format'; flush(6) - - call mesh_marc_get_tableStyles(initialcondTableStyle,hypoelasticTableStyle,FILEUNIT) - if (myDebug) write(6,'(a)') ' Got table styles'; flush(6) - - if (fileFormatVersion > 12) then - Marc_matNumber = mesh_marc_get_matNumber(FILEUNIT,hypoelasticTableStyle) - if (myDebug) write(6,'(a)') ' Got hypoleastic material number'; flush(6) - endif - - call mesh_marc_count_nodesAndElements(mesh_nNodes, mesh_nElems, FILEUNIT) - if (myDebug) write(6,'(a)') ' Counted nodes/elements'; flush(6) - - call mesh_marc_count_elementSets(mesh_NelemSets,mesh_maxNelemInSet,FILEUNIT) - if (myDebug) write(6,'(a)') ' Counted element sets'; flush(6) + allocate(calcMode(elem%nIPs,nElems),source=.false.) ! pretend to have collected what first call is asking (F = I) + calcMode(ip,mesh_FEasCP('elem',el)) = .true. ! first ip,el needs to be already pingponged to "calc" - allocate(mesh_nameElemSet(mesh_NelemSets)); mesh_nameElemSet = 'n/a' - allocate(mesh_mapElemSet(1+mesh_maxNelemInSet,mesh_NelemSets),source=0) - call mesh_marc_map_elementSets(mesh_nameElemSet,mesh_mapElemSet,FILEUNIT) - if (myDebug) write(6,'(a)') ' Mapped element sets'; flush(6) + + allocate(mesh_ipCoordinates(3,elem%nIPs,nElems),source=0.0_pReal) ! deprecated + + allocate(cellNodeDefinition(elem%nNodes-1)) + allocate(connectivity_cell(elem%NcellNodesPerCell,elem%nIPs,nElems)) + call buildCells(connectivity_cell,cellNodeDefinition,& + elem,connectivity_elem) + allocate(node0_cell(3,maxval(connectivity_cell))) + call buildCellNodes(node0_cell,& + cellNodeDefinition,node0_elem) + allocate(ip_reshaped(3,elem%nIPs*nElems),source=0.0_pReal) + call buildIPcoordinates(ip_reshaped,reshape(connectivity_cell,[elem%NcellNodesPerCell,& + elem%nIPs*nElems]),node0_cell) + + call discretization_init(microstructureAt,homogenizationAt,& + ip_reshaped,& + node0_elem) + + call writeGeometry(0,connectivity_elem,& + reshape(connectivity_cell,[elem%NcellNodesPerCell,elem%nIPs*nElems]),& + node0_cell,ip_reshaped) + + call geometry_plastic_nonlocal_setIPvolume(IPvolume(elem,node0_cell,connectivity_cell)) + unscaledNormals = IPareaNormal(elem,nElems,connectivity_cell,node0_cell) + call geometry_plastic_nonlocal_setIParea(norm2(unscaledNormals,1)) + call geometry_plastic_nonlocal_setIPareaNormal(unscaledNormals/spread(norm2(unscaledNormals,1),1,3)) + call geometry_plastic_nonlocal_results - allocate (mesh_mapFEtoCPelem(2,mesh_nElems), source = 0) - call mesh_marc_map_elements(hypoelasticTableStyle,mesh_nameElemSet,mesh_mapElemSet,mesh_nElems,fileFormatVersion,FILEUNIT) - if (myDebug) write(6,'(a)') ' Mapped elements'; flush(6) - - allocate (mesh_mapFEtoCPnode(2,mesh_Nnodes),source=0) - call mesh_marc_map_nodes(mesh_Nnodes,FILEUNIT) !ToDo: don't work on global variables - if (myDebug) write(6,'(a)') ' Mapped nodes'; flush(6) - - mesh_node0 = mesh_marc_build_nodes(mesh_Nnodes,FILEUNIT) - mesh_node = mesh_node0 - if (myDebug) write(6,'(a)') ' Built nodes'; flush(6) - - elemType = mesh_marc_getElemType(mesh_nElems,FILEUNIT) - if (myDebug) write(6,'(a)') ' Counted CP sizes'; flush(6) - - call theMesh%init('mesh',elemType,mesh_node0) - call theMesh%setNelems(mesh_nElems) - - allocate(mesh_element(4+theMesh%elem%nNodes,theMesh%nElems), source=0) - mesh_element(1,:) = -1 ! DEPRECATED - mesh_element(2,:) = elemType ! DEPRECATED - - call mesh_marc_buildElements(mesh_nElems,initialcondTableStyle,FILEUNIT) - if (myDebug) write(6,'(a)') ' Built elements'; flush(6) - close (FILEUNIT) - - call mesh_build_cellconnectivity - if (myDebug) write(6,'(a)') ' Built cell connectivity'; flush(6) - mesh_cellnode = mesh_build_cellnodes() - if (myDebug) write(6,'(a)') ' Built cell nodes'; flush(6) - - allocate(mesh_ipCoordinates(3,theMesh%elem%nIPs,theMesh%nElems),source=0.0_pReal) - call mesh_build_ipCoordinates - if (myDebug) write(6,'(a)') ' Built IP coordinates'; flush(6) - call mesh_build_ipVolumes - if (myDebug) write(6,'(a)') ' Built IP volumes'; flush(6) - call mesh_build_ipAreas - if (myDebug) write(6,'(a)') ' Built IP areas'; flush(6) - - call IP_neighborhood2 - if (myDebug) write(6,'(a)') ' Built IP neighborhood'; flush(6) - - if (usePingPong .and. (mesh_Nelems /= theMesh%nElems)) & - call IO_error(600) ! ping-pong must be disabled when having non-DAMASK elements - if (debug_e < 1 .or. debug_e > theMesh%nElems) & - call IO_error(602,ext_msg='element') ! selected element does not exist - if (debug_i < 1 .or. debug_i > theMesh%elem%nIPs) & - call IO_error(602,ext_msg='IP') ! selected element does not have requested IP - - FEsolving_execElem = [ 1,theMesh%nElems ] ! parallel loop bounds set to comprise all DAMASK elements - allocate(FEsolving_execIP(2,theMesh%nElems), source=1) ! parallel loop bounds set to comprise from first IP... - FEsolving_execIP(2,:) = theMesh%elem%nIPs - - allocate(calcMode(theMesh%elem%nIPs,theMesh%nElems)) - calcMode = .false. ! pretend to have collected what first call is asking (F = I) - calcMode(ip,mesh_FEasCP('elem',el)) = .true. ! first ip,el needs to be already pingponged to "calc" - - call discretization_init(mesh_element(3,:),mesh_element(4,:),& - reshape(mesh_ipCoordinates,[3,theMesh%elem%nIPs*theMesh%nElems]),& - mesh_node0) - call geometry_plastic_nonlocal_setIPvolume(mesh_ipVolume) - call geometry_plastic_nonlocal_setIPneighborhood(mesh_ipNeighborhood2) - call geometry_plastic_nonlocal_setIParea(mesh_IParea) - call geometry_plastic_nonlocal_setIPareaNormal(mesh_IPareaNormal) end subroutine mesh_init +!-------------------------------------------------------------------------------------------------- +!> @brief Writes all information needed for the DADF5 geometry +!-------------------------------------------------------------------------------------------------- +subroutine writeGeometry(elemType, & + connectivity_elem,connectivity_cell, & + coordinates_nodes,coordinates_points) + + integer, intent(in) :: elemType + integer, dimension(:,:), intent(in) :: & + connectivity_elem, & + connectivity_cell + real(pReal), dimension(:,:), intent(in) :: & + coordinates_nodes, & + coordinates_points + + integer, dimension(:,:), allocatable :: & + connectivity_temp + real(pReal), dimension(:,:), allocatable :: & + coordinates_temp + +#if defined(DAMASK_HDF5) + call results_openJobFile + call HDF5_closeGroup(results_addGroup('geometry')) + + connectivity_temp = connectivity_elem + call results_writeDataset('geometry',connectivity_temp,'T_e',& + 'connectivity of the elements','-') + + connectivity_temp = connectivity_cell + call results_writeDataset('geometry',connectivity_temp,'T_c', & + 'connectivity of the cells','-') + + coordinates_temp = coordinates_nodes + call results_writeDataset('geometry',coordinates_temp,'x_n', & + 'coordinates of the nodes','m') + + coordinates_temp = coordinates_points + call results_writeDataset('geometry',coordinates_temp,'x_p', & + 'coordinates of the material points','m') + + call results_closeJobFile +#endif + +end subroutine writeGeometry + + +subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogenizationAt) + + type(tElement), intent(out) :: elem + real(pReal), dimension(:,:), allocatable, intent(out) :: & + node0_elem !< node x,y,z coordinates (initially!) + integer, dimension(:,:), allocatable, intent(out) :: & + connectivity_elem + + integer, dimension(:), allocatable, intent(out) :: & + microstructureAt, & + homogenizationAt + + integer :: & + fileFormatVersion, & + hypoelasticTableStyle, & + initialcondTableStyle, & + nNodes, & + nElems + integer, parameter :: & + FILEUNIT = 222 + integer, dimension(:), allocatable :: & + matNumber !< material numbers for hypoelastic material + character(len=pStringLen), dimension(:), allocatable :: inputFile !< file content, separated per lines + + character(len=64), dimension(:), allocatable :: & + nameElemSet + integer, dimension(:,:), allocatable :: & + mapElemSet !< list of elements in elementSet + + inputFile = IO_read_ASCII(trim(modelName)//trim(InputFileExtension)) + call inputRead_fileFormat(fileFormatVersion, & + inputFile) + call inputRead_tableStyles(initialcondTableStyle,hypoelasticTableStyle, & + inputFile) + if (fileFormatVersion > 12) & + call inputRead_matNumber(matNumber, & + hypoelasticTableStyle,inputFile) + call inputRead_NnodesAndElements(nNodes,nElems,& + inputFile) + + call IO_open_inputFile(FILEUNIT,modelName) ! ToDo: It would be better to use fileContent + + call inputRead_mapElemSets(nameElemSet,mapElemSet,& + FILEUNIT) + + allocate (mesh_mapFEtoCPelem(2,nElems), source=0) + call inputRead_mapElems(hypoelasticTableStyle,nameElemSet,mapElemSet,fileFormatVersion,matNumber,FILEUNIT) + + allocate (mesh_mapFEtoCPnode(2,Nnodes), source=0) + call inputRead_mapNodes(inputFile) + + call inputRead_elemType(elem, & + nElems,inputFile) + call inputRead_elemNodes(node0_elem, & + Nnodes,inputFile) + + connectivity_elem = inputRead_connectivityElem(nElems,elem%nNodes,inputFile) + + call inputRead_microstructureAndHomogenization(microstructureAt,homogenizationAt, & + nElems,elem%nNodes,nameElemSet,mapElemSet,& + initialcondTableStyle,FILEUNIT) + close(FILEUNIT) + +end subroutine inputRead + + !-------------------------------------------------------------------------------------------------- !> @brief Figures out version of Marc input file format !-------------------------------------------------------------------------------------------------- -integer function mesh_marc_get_fileFormat(fileUnit) +subroutine inputRead_fileFormat(fileFormat,fileContent) - integer, intent(in) :: fileUnit - + integer, intent(out) :: fileFormat + character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines + integer, allocatable, dimension(:) :: chunkPos - character(len=300) line - - - rewind(fileUnit) - do - read (fileUnit,'(A300)',END=620) line - chunkPos = IO_stringPos(line) - - if ( IO_lc(IO_stringValue(line,chunkPos,1)) == 'version') then - mesh_marc_get_fileFormat = IO_intValue(line,chunkPos,2) + integer :: l + + do l = 1, size(fileContent) + chunkPos = IO_stringPos(fileContent(l)) + if ( IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'version') then + fileFormat = IO_intValue(fileContent(l),chunkPos,2) exit endif enddo -620 end function mesh_marc_get_fileFormat +end subroutine inputRead_fileFormat !-------------------------------------------------------------------------------------------------- !> @brief Figures out table styles for initial cond and hypoelastic !-------------------------------------------------------------------------------------------------- -subroutine mesh_marc_get_tableStyles(initialcond,hypoelastic,fileUnit) +subroutine inputRead_tableStyles(initialcond,hypoelastic,fileContent) - integer, intent(out) :: initialcond, hypoelastic - integer, intent(in) :: fileUnit + integer, intent(out) :: initialcond, hypoelastic + character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines integer, allocatable, dimension(:) :: chunkPos - character(len=300) line - + integer :: l + initialcond = 0 hypoelastic = 0 - rewind(fileUnit) - do - read (fileUnit,'(A300)',END=620) line - chunkPos = IO_stringPos(line) - - if ( IO_lc(IO_stringValue(line,chunkPos,1)) == 'table' .and. chunkPos(1) > 5) then - initialcond = IO_intValue(line,chunkPos,4) - hypoelastic = IO_intValue(line,chunkPos,5) + do l = 1, size(fileContent) + chunkPos = IO_stringPos(fileContent(l)) + if ( IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'table' .and. chunkPos(1) > 5) then + initialcond = IO_intValue(fileContent(l),chunkPos,4) + hypoelastic = IO_intValue(fileContent(l),chunkPos,5) exit endif enddo -620 end subroutine mesh_marc_get_tableStyles +end subroutine inputRead_tableStyles !-------------------------------------------------------------------------------------------------- !> @brief Figures out material number of hypoelastic material !-------------------------------------------------------------------------------------------------- -function mesh_marc_get_matNumber(fileUnit,tableStyle) +subroutine inputRead_matNumber(matNumber, & + tableStyle,fileContent) - integer, intent(in) :: fileUnit, tableStyle - integer, dimension(:), allocatable :: mesh_marc_get_matNumber + integer, allocatable, dimension(:), intent(out) :: matNumber + integer, intent(in) :: tableStyle + character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines + integer, allocatable, dimension(:) :: chunkPos - integer :: i, j, data_blocks - character(len=300) :: line + integer :: i, j, data_blocks, l - data_blocks = 1 - rewind(fileUnit) - do - read (fileUnit,'(A300)',END=620) line - chunkPos = IO_stringPos(line) - - if ( IO_lc(IO_stringValue(line,chunkPos,1)) == 'hypoelastic') then - read (fileUnit,'(A300)',END=620) line - if (len(trim(line))/=0) then - chunkPos = IO_stringPos(line) - data_blocks = IO_intValue(line,chunkPos,1) - endif - allocate(mesh_marc_get_matNumber(data_blocks), source = 0) - do i=1,data_blocks ! read all data blocks - read (fileUnit,'(A300)',END=620) line - chunkPos = IO_stringPos(line) - mesh_marc_get_matNumber(i) = IO_intValue(line,chunkPos,1) - do j=1,2 + tableStyle ! read 2 or 3 remaining lines of data block - read (fileUnit,'(A300)') line - enddo - enddo - exit - endif + do l = 1, size(fileContent) + chunkPos = IO_stringPos(fileContent(l)) + if ( IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'hypoelastic') then + if (len(trim(fileContent(l+1)))/=0) then + chunkPos = IO_stringPos(fileContent(l+1)) + data_blocks = IO_intValue(fileContent(l+1),chunkPos,1) + else + data_blocks = 1 + endif + allocate(matNumber(data_blocks), source = 0) + do i = 0, data_blocks - 1 + j = i*(2+tableStyle) + 1 + chunkPos = IO_stringPos(fileContent(l+1+j)) + matNumber(i+1) = IO_intValue(fileContent(l+1+j),chunkPos,1) + enddo + exit + endif enddo -620 end function mesh_marc_get_matNumber +end subroutine inputRead_matNumber !-------------------------------------------------------------------------------------------------- !> @brief Count overall number of nodes and elements !-------------------------------------------------------------------------------------------------- -subroutine mesh_marc_count_nodesAndElements(nNodes, nElems, fileUnit) - - integer, intent(in) :: fileUnit - integer, intent(out) :: nNodes, nElems +subroutine inputRead_NnodesAndElements(nNodes,nElems,& + fileContent) + integer, intent(out) :: nNodes, nElems + character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines + integer, allocatable, dimension(:) :: chunkPos - character(len=300) :: line + integer :: l nNodes = 0 nElems = 0 - rewind(fileUnit) - do - read (fileUnit,'(A300)',END=620) line - chunkPos = IO_stringPos(line) - - if ( IO_lc(IO_StringValue(line,chunkPos,1)) == 'sizing') & - nElems = IO_IntValue (line,chunkPos,3) - if ( IO_lc(IO_StringValue(line,chunkPos,1)) == 'coordinates') then - read (fileUnit,'(A300)') line - chunkPos = IO_stringPos(line) - nNodes = IO_IntValue (line,chunkPos,2) - exit ! assumes that "coordinates" comes later in file + do l = 1, size(fileContent) + chunkPos = IO_stringPos(fileContent(l)) + if (IO_lc(IO_StringValue(fileContent(l),chunkPos,1)) == 'sizing') then + nElems = IO_IntValue (fileContent(l),chunkPos,3) + elseif (IO_lc(IO_StringValue(fileContent(l),chunkPos,1)) == 'coordinates') then + chunkPos = IO_stringPos(fileContent(l+1)) + nNodes = IO_IntValue (fileContent(l+1),chunkPos,2) endif enddo -620 end subroutine mesh_marc_count_nodesAndElements +end subroutine inputRead_NnodesAndElements !-------------------------------------------------------------------------------------------------- !> @brief Count overall number of element sets in mesh. !-------------------------------------------------------------------------------------------------- -subroutine mesh_marc_count_elementSets(nElemSets,maxNelemInSet,fileUnit) +subroutine inputRead_NelemSets(nElemSets,maxNelemInSet,& + fileUnit) integer, intent(out) :: nElemSets, maxNelemInSet integer, intent(in) :: fileUnit @@ -405,22 +380,26 @@ subroutine mesh_marc_count_elementSets(nElemSets,maxNelemInSet,fileUnit) endif enddo -620 end subroutine mesh_marc_count_elementSets +620 end subroutine inputRead_NelemSets !-------------------------------------------------------------------------------------------------- !> @brief map element sets !-------------------------------------------------------------------------------------------------- -subroutine mesh_marc_map_elementSets(nameElemSet,mapElemSet,fileUnit) +subroutine inputRead_mapElemSets(nameElemSet,mapElemSet,fileUnit) - character(len=64), dimension(:), intent(out) :: nameElemSet - integer, dimension(:,:), intent(out) :: mapElemSet - integer, intent(in) :: fileUnit + character(len=64), dimension(:), allocatable, intent(out) :: nameElemSet + integer, dimension(:,:), allocatable, intent(out) :: mapElemSet + integer, intent(in) :: fileUnit integer, allocatable, dimension(:) :: chunkPos character(len=300) :: line - integer :: elemSet - + integer :: elemSet, NelemSets, maxNelemInSet + + + call inputRead_NelemSets(NelemSets,maxNelemInSet,fileUnit) + allocate(nameElemSet(NelemSets)); nameElemSet = 'n/a' + allocate(mapElemSet(1+maxNelemInSet,NelemSets),source=0) elemSet = 0 rewind(fileUnit) @@ -435,179 +414,172 @@ subroutine mesh_marc_map_elementSets(nameElemSet,mapElemSet,fileUnit) endif enddo -620 end subroutine mesh_marc_map_elementSets +620 end subroutine inputRead_mapElemSets !-------------------------------------------------------------------------------------------------- !> @brief Maps elements from FE ID to internal (consecutive) representation. !-------------------------------------------------------------------------------------------------- -subroutine mesh_marc_map_elements(tableStyle,nameElemSet,mapElemSet,nElems,fileFormatVersion,fileUnit) +subroutine inputRead_mapElems(tableStyle,nameElemSet,mapElemSet,fileFormatVersion,matNumber,fileUnit) - integer, intent(in) :: fileUnit,tableStyle,nElems,fileFormatVersion - character(len=64), intent(in), dimension(:) :: nameElemSet - integer, dimension(:,:), intent(in) :: & - mapElemSet - - integer, allocatable, dimension(:) :: chunkPos - character(len=300) :: line, & - tmp - - integer, dimension (1+nElems) :: contInts - integer :: i,cpElem - - cpElem = 0 - contInts = 0 - rewind(fileUnit) - do - read (fileUnit,'(A300)',END=620) line - chunkPos = IO_stringPos(line) - if (fileFormatVersion < 13) then ! Marc 2016 or earlier - if( IO_lc(IO_stringValue(line,chunkPos,1)) == 'hypoelastic' ) then - do i=1,3+TableStyle ! skip three (or four if new table style!) lines - read (fileUnit,'(A300)') line - enddo - contInts = IO_continuousIntValues(fileUnit,nElems,nameElemSet,& - mapElemSet,size(nameElemSet)) - exit - endif - else ! Marc2017 and later - if ( IO_lc(IO_stringValue(line,chunkPos,1)) == 'connectivity') then - read (fileUnit,'(A300)',END=620) line - chunkPos = IO_stringPos(line) - if(any(Marc_matNumber==IO_intValue(line,chunkPos,6))) then - do - read (fileUnit,'(A300)',END=620) line - chunkPos = IO_stringPos(line) - tmp = IO_lc(IO_stringValue(line,chunkPos,1)) - if (verify(trim(tmp),"0123456789")/=0) then ! found keyword - exit - else - contInts(1) = contInts(1) + 1 - read (tmp,*) contInts(contInts(1)+1) - endif - enddo - endif - endif - endif - enddo + integer, intent(in) :: fileUnit,tableStyle,fileFormatVersion + integer, dimension(:), intent(in) :: matNumber + character(len=64), intent(in), dimension(:) :: nameElemSet + integer, dimension(:,:), intent(in) :: & + mapElemSet + + integer, allocatable, dimension(:) :: chunkPos + character(len=300) :: line, & + tmp + + integer, dimension(:), allocatable :: contInts + integer :: i,cpElem + + allocate(contInts(size(mesh_mapFEtoCPelem,2)+1)) + + cpElem = 0 + contInts = 0 + rewind(fileUnit) + do + read (fileUnit,'(A300)',END=620) line + chunkPos = IO_stringPos(line) + Marc2016andEarlier: if (fileFormatVersion < 13) then + if( IO_lc(IO_stringValue(line,chunkPos,1)) == 'hypoelastic' ) then + skipLines: do i=1,3+TableStyle + read (fileUnit,'(A300)') line + enddo skipLines + contInts = IO_continuousIntValues(fileUnit,size(mesh_mapFEtoCPelem,2),nameElemSet,& + mapElemSet,size(nameElemSet)) + exit + endif + else Marc2016andEarlier + if ( IO_lc(IO_stringValue(line,chunkPos,1)) == 'connectivity') then + read (fileUnit,'(A300)',END=620) line + chunkPos = IO_stringPos(line) + if(any(matNumber==IO_intValue(line,chunkPos,6))) then + do + read (fileUnit,'(A300)',END=620) line + chunkPos = IO_stringPos(line) + tmp = IO_lc(IO_stringValue(line,chunkPos,1)) + if (verify(trim(tmp),"0123456789")/=0) then ! found keyword + exit + else + contInts(1) = contInts(1) + 1 + read (tmp,*) contInts(contInts(1)+1) + endif + enddo + endif + endif + endif Marc2016andEarlier + enddo 620 do i = 1,contInts(1) - cpElem = cpElem+1 - mesh_mapFEtoCPelem(1,cpElem) = contInts(1+i) - mesh_mapFEtoCPelem(2,cpElem) = cpElem - enddo + cpElem = cpElem+1 + mesh_mapFEtoCPelem(1,cpElem) = contInts(1+i) + mesh_mapFEtoCPelem(2,cpElem) = cpElem + enddo call math_sort(mesh_mapFEtoCPelem) -end subroutine mesh_marc_map_elements +end subroutine inputRead_mapElems !-------------------------------------------------------------------------------------------------- !> @brief Maps node from FE ID to internal (consecutive) representation. !-------------------------------------------------------------------------------------------------- -subroutine mesh_marc_map_nodes(nNodes,fileUnit) +subroutine inputRead_mapNodes(fileContent) - integer, intent(in) :: fileUnit, nNodes + character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines integer, allocatable, dimension(:) :: chunkPos - character(len=300) line + integer :: i, l - integer, dimension (nNodes) :: node_count - integer :: i - - node_count = 0 - - rewind(fileUnit) - do - read (fileUnit,'(A300)',END=620) line - chunkPos = IO_stringPos(line) - if( IO_lc(IO_stringValue(line,chunkPos,1)) == 'coordinates' ) then - read (fileUnit,'(A300)') line ! skip crap line - do i = 1,nNodes - read (fileUnit,'(A300)') line - mesh_mapFEtoCPnode(1,i) = IO_fixedIntValue (line,[0,10],1) + do l = 1, size(fileContent) + chunkPos = IO_stringPos(fileContent(l)) + if( IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'coordinates' ) then + do i = 1,size(mesh_mapFEtoCPnode,2) + mesh_mapFEtoCPnode(1,i) = IO_fixedIntValue (fileContent(l+1+i),[0,10],1) mesh_mapFEtoCPnode(2,i) = i enddo exit endif enddo -620 call math_sort(mesh_mapFEtoCPnode) + call math_sort(mesh_mapFEtoCPnode) -end subroutine mesh_marc_map_nodes +end subroutine inputRead_mapNodes !-------------------------------------------------------------------------------------------------- !> @brief store x,y,z coordinates of all nodes in mesh. !-------------------------------------------------------------------------------------------------- -function mesh_marc_build_nodes(nNode,fileUnit) result(nodes) +subroutine inputRead_elemNodes(nodes, & + nNode,fileContent) - integer, intent(in) :: nNode,fileUnit - real(pReal), dimension(3,nNode) :: nodes + real(pReal), allocatable, dimension(:,:), intent(out) :: nodes + integer, intent(in) :: nNode + character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines + integer, dimension(5), parameter :: node_ends = [0,10,30,50,70] integer, allocatable, dimension(:) :: chunkPos - character(len=300) :: line - integer :: i,j,m + integer :: i,j,m,l - rewind(fileUnit) - do - read (fileUnit,'(A300)',END=620) line - chunkPos = IO_stringPos(line) - if( IO_lc(IO_stringValue(line,chunkPos,1)) == 'coordinates' ) then - read (fileUnit,'(A300)') line ! skip crap line + allocate(nodes(3,nNode)) + + do l = 1, size(fileContent) + chunkPos = IO_stringPos(fileContent(l)) + if( IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'coordinates' ) then do i=1,nNode - read (fileUnit,'(A300)') line - m = mesh_FEasCP('node',IO_fixedIntValue(line,node_ends,1)) + m = mesh_FEasCP('node',IO_fixedIntValue(fileContent(l+1+i),node_ends,1)) do j = 1,3 - nodes(j,m) = mesh_unitlength * IO_fixedNoEFloatValue(line,node_ends,j+1) + nodes(j,m) = mesh_unitlength * IO_fixedNoEFloatValue(fileContent(l+1+i),node_ends,j+1) enddo enddo exit endif enddo -620 end function mesh_marc_build_nodes +end subroutine inputRead_elemNodes !-------------------------------------------------------------------------------------------------- !> @brief Gets element type (and checks if the whole mesh comprises of only one type) !-------------------------------------------------------------------------------------------------- -integer function mesh_marc_getElemType(nElem,fileUnit) +subroutine inputRead_elemType(elem, & + nElem,fileContent) - integer, intent(in) :: & - nElem, & - fileUnit + type(tElement), intent(out) :: elem + integer, intent(in) :: nElem + character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines - type(tElement) :: tempEl integer, allocatable, dimension(:) :: chunkPos - character(len=300) :: line - integer :: i,t + integer :: i,j,t,l,remainingChunks t = -1 - - rewind(fileUnit) - do - read (fileUnit,'(A300)',END=620) line - chunkPos = IO_stringPos(line) - if( IO_lc(IO_stringValue(line,chunkPos,1)) == 'connectivity' ) then - read (fileUnit,'(A300)') line ! Garbage line + do l = 1, size(fileContent) + chunkPos = IO_stringPos(fileContent(l)) + if( IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'connectivity' ) then + j = 0 do i=1,nElem ! read all elements - read (fileUnit,'(A300)') line - chunkPos = IO_stringPos(line) + chunkPos = IO_stringPos(fileContent(l+1+i+j)) if (t == -1) then - t = mapElemtype(IO_stringValue(line,chunkPos,2)) - call tempEl%init(t) - mesh_marc_getElemType = t + t = mapElemtype(IO_stringValue(fileContent(l+1+i+j),chunkPos,2)) + call elem%init(t) else - if (t /= mapElemtype(IO_stringValue(line,chunkPos,2))) call IO_error(191,el=t,ip=i) + if (t /= mapElemtype(IO_stringValue(fileContent(l+1+i+j),chunkPos,2))) call IO_error(191,el=t,ip=i) endif - call IO_skipChunks(fileUnit,tempEl%nNodes-(chunkPos(1)-2)) + remainingChunks = elem%nNodes - (chunkPos(1) - 2) + do while(remainingChunks > 0) + j = j + 1 + chunkPos = IO_stringPos(fileContent(l+1+i+j)) + remainingChunks = remainingChunks - chunkPos(1) + enddo enddo exit endif enddo -contains + contains !-------------------------------------------------------------------------------------------------- !> @brief mapping of Marc element types to internal representation @@ -650,49 +622,49 @@ contains call IO_error(error_ID=190,ext_msg=IO_lc(what)) end select -end function mapElemtype + end function mapElemtype -620 end function mesh_marc_getElemType +end subroutine inputRead_elemType !-------------------------------------------------------------------------------------------------- -!> @brief Stores node IDs and homogenization and microstructure ID +!> @brief Stores node IDs !-------------------------------------------------------------------------------------------------- -subroutine mesh_marc_buildElements(nElem,initialcondTableStyle,fileUnit) +function inputRead_connectivityElem(nElem,nNodes,fileContent) integer, intent(in) :: & nElem, & - initialcondTableStyle, & - fileUnit + nNodes !< number of nodes per element + character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines + + integer, dimension(nNodes,nElem) :: & + inputRead_connectivityElem integer, allocatable, dimension(:) :: chunkPos - character(len=300) line integer, dimension(1+nElem) :: contInts - integer :: i,j,t,sv,myVal,e,nNodesAlreadyRead + integer :: i,k,j,t,e,l,nNodesAlreadyRead - rewind(fileUnit) - do - read (fileUnit,'(A300)',END=620) line - chunkPos = IO_stringPos(line) - if( IO_lc(IO_stringValue(line,chunkPos,1)) == 'connectivity' ) then - read (fileUnit,'(A300)',END=620) line ! garbage line + do l = 1, size(fileContent) + chunkPos = IO_stringPos(fileContent(l)) + if( IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'connectivity' ) then + j = 0 do i = 1,nElem - read (fileUnit,'(A300)',END=620) line - chunkPos = IO_stringPos(line) - e = mesh_FEasCP('elem',IO_intValue(line,chunkPos,1)) + chunkPos = IO_stringPos(fileContent(l+1+i+j)) + e = mesh_FEasCP('elem',IO_intValue(fileContent(l+1+i+j),chunkPos,1)) if (e /= 0) then ! disregard non CP elems - nNodesAlreadyRead = 0 - do j = 1,chunkPos(1)-2 - mesh_element(4+j,e) = mesh_FEasCP('node',IO_IntValue(line,chunkPos,j+2)) ! CP ids of nodes + do k = 1,chunkPos(1)-2 + inputRead_connectivityElem(k,e) = & + mesh_FEasCP('node',IO_IntValue(fileContent(l+1+i+j),chunkPos,k+2)) enddo nNodesAlreadyRead = chunkPos(1) - 2 - do while(nNodesAlreadyRead < theMesh%elem%nNodes) ! read on if not all nodes in one line - read (fileUnit,'(A300)',END=620) line - chunkPos = IO_stringPos(line) - do j = 1,chunkPos(1) - mesh_element(4+nNodesAlreadyRead+j,e) = mesh_FEasCP('node',IO_IntValue(line,chunkPos,j)) ! CP ids of nodes + do while(nNodesAlreadyRead < nNodes) ! read on if not all nodes in one line + j = j + 1 + chunkPos = IO_stringPos(fileContent(l+1+i+j)) + do k = 1,chunkPos(1) + inputRead_connectivityElem(nNodesAlreadyRead+k,e) = & + mesh_FEasCP('node',IO_IntValue(fileContent(l+1+i+j),chunkPos,k)) enddo nNodesAlreadyRead = nNodesAlreadyRead + chunkPos(1) enddo @@ -701,18 +673,40 @@ subroutine mesh_marc_buildElements(nElem,initialcondTableStyle,fileUnit) exit endif enddo -620 rewind(fileUnit) ! just in case "initial state" appears before "connectivity" -#if defined(DAMASK_HDF5) - call results_openJobFile - call HDF5_closeGroup(results_addGroup('geometry')) - call results_writeDataset('geometry',mesh_element(5:,:),'C',& - 'connectivity of the elements','-') - call results_closeJobFile -#endif - - call buildCells(theMesh,theMesh%elem,mesh_element(5:,:)) +end function inputRead_connectivityElem + + +!-------------------------------------------------------------------------------------------------- +!> @brief Stores homogenization and microstructure ID +!-------------------------------------------------------------------------------------------------- +subroutine inputRead_microstructureAndHomogenization(microstructureAt,homogenizationAt, & + nElem,nNodes,nameElemSet,mapElemSet,initialcondTableStyle,fileUnit) + + integer, dimension(:), allocatable, intent(out) :: & + microstructureAt, & + homogenizationAt + integer, intent(in) :: & + nElem, & + nNodes, & !< number of nodes per element + initialcondTableStyle, & + fileUnit + character(len=64), dimension(:), intent(in) :: & + nameElemSet + integer, dimension(:,:), intent(in) :: & + mapElemSet !< list of elements in elementSet + + integer, allocatable, dimension(:) :: chunkPos + character(len=300) line + integer, dimension(1+nElem) :: contInts + integer :: i,j,t,sv,myVal,e,nNodesAlreadyRead + + + allocate(microstructureAt(nElem),source=0) + allocate(homogenizationAt(nElem),source=0) + + rewind(fileUnit) read (fileUnit,'(A300)',END=630) line do chunkPos = IO_stringPos(line) @@ -732,10 +726,11 @@ subroutine mesh_marc_buildElements(nElem,initialcondTableStyle,fileUnit) read (fileUnit,'(A300)',END=630) line ! read extra line endif contInts = IO_continuousIntValues& ! get affected elements - (fileUnit,theMesh%nElems,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets) + (fileUnit,nElem,nameElemSet,mapElemSet,size(nameElemSet)) do i = 1,contInts(1) e = mesh_FEasCP('elem',contInts(1+i)) - mesh_element(1+sv,e) = myVal + if (sv == 2) microstructureAt(e) = myVal + if (sv == 3) homogenizationAt(e) = myVal enddo if (initialcondTableStyle == 0) read (fileUnit,'(A300)',END=630) line ! ignore IP range for old table style read (fileUnit,'(A300)',END=630) line @@ -747,26 +742,30 @@ subroutine mesh_marc_buildElements(nElem,initialcondTableStyle,fileUnit) endif enddo -630 end subroutine mesh_marc_buildElements +630 end subroutine inputRead_microstructureAndHomogenization -subroutine buildCells(thisMesh,elem,connectivity_elem) +!-------------------------------------------------------------------------------------------------- +!> @brief Calculates cell node coordinates from element node coordinates +!-------------------------------------------------------------------------------------------------- +subroutine buildCells(connectivity_cell,cellNodeDefinition, & + elem,connectivity_elem) - class(tMesh) :: thisMesh - type(tElement) :: elem - integer,dimension(:,:), intent(in) :: connectivity_elem - integer,dimension(:,:), allocatable :: parentsAndWeights,candidates_global - integer,dimension(:), allocatable :: candidates_local - integer,dimension(:,:,:), allocatable :: connectivity_cell - integer,dimension(:,:), allocatable :: connectivity_cell_reshape - real(pReal), dimension(:,:), allocatable :: nodes_new,nodes + type(tCellNodeDefinition), dimension(:), intent(out) :: cellNodeDefinition ! definition of cell nodes for increasing number of parents + integer, dimension(:,:,:),intent(out) :: connectivity_cell + + type(tElement), intent(in) :: elem ! element definition + integer, dimension(:,:), intent(in) :: connectivity_elem ! connectivity of the elements + + integer,dimension(:), allocatable :: candidates_local + integer,dimension(:,:), allocatable :: parentsAndWeights,candidates_global + integer :: e, n, c, p, s,i,m,j,nParentNodes,nCellNode,Nelem,candidateID - Nelem = thisMesh%Nelems + Nelem = size(connectivity_elem,2) !--------------------------------------------------------------------------------------------------- ! initialize global connectivity to negative local connectivity - allocate(connectivity_cell(elem%NcellNodesPerCell,elem%nIPs,Nelem)) connectivity_cell = -spread(elem%cell,3,Nelem) ! local cell node ID !--------------------------------------------------------------------------------------------------- @@ -781,8 +780,8 @@ subroutine buildCells(thisMesh,elem,connectivity_elem) endif realNode enddo enddo - nCellNode = thisMesh%nNodes - + + nCellNode = maxval(connectivity_elem) !--------------------------------------------------------------------------------------------------- ! set connectivity of cell nodes that are defined by 2,...,nNodes real nodes @@ -843,53 +842,37 @@ subroutine buildCells(thisMesh,elem,connectivity_elem) enddo i = uniqueRows(candidates_global(1:2*nParentNodes,:)) - - - ! calculate coordinates of cell nodes and insert their ID into the cell conectivity - nodes_new = reshape([(0.0_pReal,j = 1, 3*i)], [3,i]) + allocate(cellNodeDefinition(nParentNodes-1)%parents(i,nParentNodes)) + allocate(cellNodeDefinition(nParentNodes-1)%weights(i,nParentNodes)) i = 1 n = 1 do while(n <= size(candidates_local)*Nelem) - j=0 - parentsAndWeights(:,1) = candidates_global(1:nParentNodes,n+j) - parentsAndWeights(:,2) = candidates_global(nParentNodes+1:nParentNodes*2,n+j) - e = candidates_global(nParentNodes*2+1,n+j) - c = candidates_global(nParentNodes*2+2,n+j) - do m = 1, nParentNodes - nodes_new(:,i) = nodes_new(:,i) & - + thisMesh%node_0(:,parentsAndWeights(m,1)) * real(parentsAndWeights(m,2),pReal) - enddo - nodes_new(:,i) = nodes_new(:,i)/real(sum(parentsAndWeights(:,2)),pReal) + j=0 + parentsAndWeights(:,1) = candidates_global(1:nParentNodes,n+j) + parentsAndWeights(:,2) = candidates_global(nParentNodes+1:nParentNodes*2,n+j) - do while (n+j<= size(candidates_local)*Nelem) - if (any(candidates_global(1:2*nParentNodes,n+j)/=candidates_global(1:2*nParentNodes,n))) exit - where (connectivity_cell(:,:,candidates_global(nParentNodes*2+1,n+j)) == -candidates_global(nParentNodes*2+2,n+j)) ! still locally defined - connectivity_cell(:,:,candidates_global(nParentNodes*2+1,n+j)) = nCellNode + i ! gets current new cell node id - end where + e = candidates_global(nParentNodes*2+1,n+j) + c = candidates_global(nParentNodes*2+2,n+j) + + do while (n+j<= size(candidates_local)*Nelem) + if (any(candidates_global(1:2*nParentNodes,n+j)/=candidates_global(1:2*nParentNodes,n))) exit + where (connectivity_cell(:,:,candidates_global(nParentNodes*2+1,n+j)) == -candidates_global(nParentNodes*2+2,n+j)) ! still locally defined + connectivity_cell(:,:,candidates_global(nParentNodes*2+1,n+j)) = nCellNode + 1 ! gets current new cell node id + end where - j = j + 1 - enddo - i=i+1 - n = n+j - + j = j+1 + enddo + nCellNode = nCellNode + 1 + cellNodeDefinition(nParentNodes-1)%parents(i,:) = parentsAndWeights(:,1) + cellNodeDefinition(nParentNodes-1)%weights(i,:) = parentsAndWeights(:,2) + i = i + 1 + n = n+j enddo - nCellNode = nCellNode + i - if (i/=0) nodes = reshape([nodes,nodes_new],[3,nCellNode]) + enddo - thisMesh%node_0 = nodes - mesh_cell2 = connectivity_cell - -#if defined(DAMASK_HDF5) - connectivity_cell_reshape = reshape(connectivity_cell,[elem%NcellNodesPerCell,elem%nIPs*Nelem]) - call results_openJobFile - call results_writeDataset('geometry',connectivity_cell_reshape,'c',& - 'connectivity of the cells','-') - call results_closeJobFile -#endif contains - !------------------------------------------------------------------------------------------------ !> @brief count unique rows (same rows need to be stored consecutively) !------------------------------------------------------------------------------------------------ @@ -920,340 +903,166 @@ end subroutine buildCells !-------------------------------------------------------------------------------------------------- -!> @brief Split CP elements into cells. -!> @details Build a mapping between cells and the corresponding cell nodes ('mesh_cell'). -!> Cell nodes that are also matching nodes are unique in the list of cell nodes, -!> all others (currently) might be stored more than once. +!> @brief Calculates cell node coordinates from element node coordinates !-------------------------------------------------------------------------------------------------- -subroutine mesh_build_cellconnectivity +subroutine buildCellNodes(node_cell, & + definition,node_elem) - integer, dimension(:), allocatable :: & - matchingNode2cellnode - integer, dimension(:,:), allocatable :: & - cellnodeParent - integer, dimension(theMesh%elem%Ncellnodes) :: & - localCellnode2globalCellnode - integer :: & - e,n,i, & - matchingNodeID, & - localCellnodeID + real(pReal), dimension(:,:), intent(out) :: node_cell !< cell node coordinates + type(tCellNodeDefinition), dimension(:), intent(in) :: definition !< cell node definition (weights and parents) + real(pReal), dimension(:,:), intent(in) :: node_elem !< element nodes + + integer :: i, j, k, n + + n = size(node_elem,2) + node_cell(:,1:n) = node_elem !< initial nodes coincide with element nodes - allocate(mesh_cell(FE_maxNcellnodesPerCell,theMesh%elem%nIPs,theMesh%nElems), source=0) - allocate(matchingNode2cellnode(theMesh%nNodes), source=0) - allocate(cellnodeParent(2,theMesh%elem%Ncellnodes*theMesh%nElems), source=0) - - mesh_Ncells = theMesh%nElems*theMesh%elem%nIPs -!-------------------------------------------------------------------------------------------------- -! Count cell nodes (including duplicates) and generate cell connectivity list - mesh_Ncellnodes = 0 - - do e = 1,theMesh%nElems - localCellnode2globalCellnode = 0 - do i = 1,theMesh%elem%nIPs - do n = 1,theMesh%elem%NcellnodesPerCell - localCellnodeID = theMesh%elem%cell(n,i) - if (localCellnodeID <= FE_NmatchingNodes(theMesh%elem%geomType)) then ! this cell node is a matching node - matchingNodeID = mesh_element(4+localCellnodeID,e) - if (matchingNode2cellnode(matchingNodeID) == 0) then ! if this matching node does not yet exist in the glbal cell node list ... - mesh_Ncellnodes = mesh_Ncellnodes + 1 ! ... count it as cell node ... - matchingNode2cellnode(matchingNodeID) = mesh_Ncellnodes ! ... and remember its global ID - cellnodeParent(1,mesh_Ncellnodes) = e ! ... and where it belongs to - cellnodeParent(2,mesh_Ncellnodes) = localCellnodeID - endif - mesh_cell(n,i,e) = matchingNode2cellnode(matchingNodeID) - else ! this cell node is no matching node - if (localCellnode2globalCellnode(localCellnodeID) == 0) then ! if this local cell node does not yet exist in the global cell node list ... - mesh_Ncellnodes = mesh_Ncellnodes + 1 ! ... count it as cell node ... - localCellnode2globalCellnode(localCellnodeID) = mesh_Ncellnodes ! ... and remember its global ID ... - cellnodeParent(1,mesh_Ncellnodes) = e ! ... and it belongs to - cellnodeParent(2,mesh_Ncellnodes) = localCellnodeID - endif - mesh_cell(n,i,e) = localCellnode2globalCellnode(localCellnodeID) - endif - enddo - enddo - enddo - - allocate(mesh_cellnodeParent(2,mesh_Ncellnodes)) - allocate(mesh_cellnode(3,mesh_Ncellnodes)) - - forall(n = 1:mesh_Ncellnodes) - mesh_cellnodeParent(1,n) = cellnodeParent(1,n) - mesh_cellnodeParent(2,n) = cellnodeParent(2,n) - endforall - -end subroutine mesh_build_cellconnectivity - - -!-------------------------------------------------------------------------------------------------- -!> @brief Calculate position of cellnodes from the given position of nodes -!> Build list of cellnodes' coordinates. -!> Cellnode coordinates are calculated from a weighted sum of node coordinates. -!-------------------------------------------------------------------------------------------------- -function mesh_build_cellnodes() - - - real(pReal), dimension(3,mesh_Ncellnodes) :: mesh_build_cellnodes - - integer :: & - e,n,m, & - localCellnodeID - real(pReal), dimension(3) :: & - myCoords - - mesh_build_cellnodes = 0.0_pReal -!$OMP PARALLEL DO PRIVATE(e,localCellnodeID,myCoords) - do n = 1,mesh_Ncellnodes ! loop over cell nodes - e = mesh_cellnodeParent(1,n) - localCellnodeID = mesh_cellnodeParent(2,n) - myCoords = 0.0_pReal - do m = 1,theMesh%elem%nNodes - myCoords = myCoords + mesh_node(1:3,mesh_element(4+m,e)) & - * theMesh%elem%cellNodeParentNodeWeights(m,localCellnodeID) - enddo - mesh_build_cellnodes(1:3,n) = myCoords / sum(theMesh%elem%cellNodeParentNodeWeights(:,localCellnodeID)) - enddo -!$OMP END PARALLEL DO - -end function mesh_build_cellnodes - - -!-------------------------------------------------------------------------------------------------- -!> @brief Calculates IP volume. Allocates global array 'mesh_ipVolume' -!> @details The IP volume is calculated differently depending on the cell type. -!> 2D cells assume an element depth of one in order to calculate the volume. -!> For the hexahedral cell we subdivide the cell into subvolumes of pyramidal -!> shape with a cell face as basis and the central ip at the tip. This subvolume is -!> calculated as an average of four tetrahedals with three corners on the cell face -!> and one corner at the central ip. -!-------------------------------------------------------------------------------------------------- -subroutine mesh_build_ipVolumes - - integer :: e,i,c,m,f,n - real(pReal), dimension(size(theMesh%elem%cellFace,1),size(theMesh%elem%cellFace,2)) :: subvolume - - allocate(mesh_ipVolume(theMesh%elem%nIPs,theMesh%nElems),source=0.0_pReal) - c = theMesh%elem%cellType - m = FE_NcellnodesPerCellface(c) - - !$OMP PARALLEL DO PRIVATE(f,n,subvolume) - do e = 1,theMesh%nElems - select case (c) - - case (1) ! 2D 3node - forall (i = 1:theMesh%elem%nIPs) & ! loop over ips=cells in this element - mesh_ipVolume(i,e) = math_areaTriangle(theMesh%node_0(1:3,mesh_cell2(1,i,e)), & - theMesh%node_0(1:3,mesh_cell2(2,i,e)), & - theMesh%node_0(1:3,mesh_cell2(3,i,e))) - - case (2) ! 2D 4node - forall (i = 1:theMesh%elem%nIPs) & ! loop over ips=cells in this element - mesh_ipVolume(i,e) = math_areaTriangle(theMesh%node_0(1:3,mesh_cell2(1,i,e)), & ! here we assume a planar shape, so division in two triangles suffices - theMesh%node_0(1:3,mesh_cell2(2,i,e)), & - theMesh%node_0(1:3,mesh_cell2(3,i,e))) & - + math_areaTriangle(theMesh%node_0(1:3,mesh_cell2(3,i,e)), & - theMesh%node_0(1:3,mesh_cell2(4,i,e)), & - theMesh%node_0(1:3,mesh_cell2(1,i,e))) - - case (3) ! 3D 4node - forall (i = 1:theMesh%elem%nIPs) & ! loop over ips=cells in this element - mesh_ipVolume(i,e) = math_volTetrahedron(theMesh%node_0(1:3,mesh_cell2(1,i,e)), & - theMesh%node_0(1:3,mesh_cell2(2,i,e)), & - theMesh%node_0(1:3,mesh_cell2(3,i,e)), & - theMesh%node_0(1:3,mesh_cell2(4,i,e))) - - case (4) ! 3D 8node - do i = 1,theMesh%elem%nIPs ! loop over ips=cells in this element - subvolume = 0.0_pReal - forall(f = 1:FE_NipNeighbors(c), n = 1:m) & - subvolume(n,f) = math_volTetrahedron(& - mesh_cellnode(1:3,mesh_cell(theMesh%elem%cellface( n ,f),i,e)), & - mesh_cellnode(1:3,mesh_cell(theMesh%elem%cellface(1+mod(n ,m),f),i,e)), & - mesh_cellnode(1:3,mesh_cell(theMesh%elem%cellface(1+mod(n+1,m),f),i,e)), & - mesh_ipCoordinates(1:3,i,e)) - mesh_ipVolume(i,e) = 0.5_pReal * sum(subvolume) ! each subvolume is based on four tetrahedrons, altough the face consists of only two triangles -> averaging factor two - enddo - - end select - enddo - !$OMP END PARALLEL DO - -end subroutine mesh_build_ipVolumes - - -subroutine IP_neighborhood2 - - integer, dimension(:,:), allocatable :: faces - integer, dimension(:), allocatable :: face - integer :: e,i,f,c,m,n,j,k,l,p, current, next,i2,e2,n2,k2 - logical :: match - allocate(faces(size(theMesh%elem%cellface,1)+3,size(theMesh%elem%cellface,2)*theMesh%elem%nIPs*theMesh%Nelems)) - - ! store cell face definitions - f = 0 - do e = 1,theMesh%nElems - do i = 1,theMesh%elem%nIPs - do n = 1, theMesh%elem%nIPneighbors - f = f + 1 - face = mesh_cell2(theMesh%elem%cellFace(:,n),i,e) - storeSorted: do j = 1, size(face) - faces(j,f) = maxval(face) - face(maxloc(face)) = -huge(1) - enddo storeSorted - faces(j:j+2,f) = [e,i,n] + do i = 1, size(cellNodeDefinition,1) + do j = 1, size(cellNodeDefinition(i)%parents,1) + n = n+1 + node_cell(:,n) = 0.0_pReal + do k = 1, size(cellNodeDefinition(i)%parents,2) + node_cell(:,n) = node_cell(:,n) & + + node_cell(:,definition(i)%parents(j,k)) * real(definition(i)%weights(j,k),pReal) enddo + node_cell(:,n) = node_cell(:,n)/real(sum(definition(i)%weights(j,:)),pReal) enddo enddo - ! sort .. - call math_sort(faces,sortDim=1) - do p = 2, size(faces,1)-2 - n = 1 - do while(n <= size(faces,2)) - j=0 - do while (n+j<= size(faces,2)) - if (faces(p-1,n+j)/=faces(p-1,n)) exit - j = j + 1 - enddo - e = n+j-1 - if (any(faces(p,n:e)/=faces(p,n))) call math_sort(faces(:,n:e),sortDim=p) - n = e+1 - enddo +end subroutine buildCellNodes + + +!-------------------------------------------------------------------------------------------------- +!> @brief Calculates IP coordinates as center of cell +!-------------------------------------------------------------------------------------------------- +subroutine buildIPcoordinates(IPcoordinates, & + connectivity_cell,node_cell) + + real(pReal), dimension(:,:), intent(out):: IPcoordinates !< cell-center/IP coordinates + integer, dimension(:,:), intent(in) :: connectivity_cell !< connectivity for each cell + real(pReal), dimension(:,:), intent(in) :: node_cell !< cell node coordinates + + integer :: i, n + + do i = 1, size(connectivity_cell,2) + IPcoordinates(:,i) = 0.0_pReal + do n = 1, size(connectivity_cell,1) + IPcoordinates(:,i) = IPcoordinates(:,i) & + + node_cell(:,connectivity_cell(n,i)) + enddo + IPcoordinates(:,i) = IPcoordinates(:,i)/real(size(connectivity_cell,1),pReal) enddo + +end subroutine buildIPcoordinates - allocate(mesh_ipNeighborhood2(3,theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems),source=0) + +!--------------------------------------------------------------------------------------------------- +!> @brief Calculates IP volume. +!> @details The IP volume is calculated differently depending on the cell type. +!> 2D cells assume an element depth of 1.0 +!--------------------------------------------------------------------------------------------------- +function IPvolume(elem,node,connectivity) - ! find IP neighbors - f = 1 - do while(f <= size(faces,2)) - e = faces(size(theMesh%elem%cellface,1)+1,f) - i = faces(size(theMesh%elem%cellface,1)+2,f) - n = faces(size(theMesh%elem%cellface,1)+3,f) - - if (f < size(faces,2)) then - match = all(faces(1:size(theMesh%elem%cellface,1),f) == faces(1:size(theMesh%elem%cellface,1),f+1)) - e2 = faces(size(theMesh%elem%cellface,1)+1,f+1) - i2 = faces(size(theMesh%elem%cellface,1)+2,f+1) - n2 = faces(size(theMesh%elem%cellface,1)+3,f+1) - else - match = .false. - endif - - if (match) then - if (e == e2) then ! same element. MD: I don't think that we need this (not even for other elements) - k = theMesh%elem%IPneighbor(n, i) - k2 = theMesh%elem%IPneighbor(n2,i2) - endif - mesh_ipNeighborhood2(1:3,n, i, e) = [e2,i2,n2] - mesh_ipNeighborhood2(1:3,n2,i2,e2) = [e, i, n] - f = f +1 - endif - f = f +1 - enddo + type(tElement), intent(in) :: elem + real(pReal), dimension(:,:), intent(in) :: node + integer, dimension(:,:,:), intent(in) :: connectivity + + real(pReal), dimension(elem%nIPs,size(connectivity,3)) :: IPvolume + real(pReal), dimension(3) :: x0,x1,x2,x3,x4,x5,x6,x7 -end subroutine IP_neighborhood2 + integer :: e,i -!-------------------------------------------------------------------------------------------------- -!> @brief Calculates IP Coordinates. Allocates global array 'mesh_ipCoordinates' -! Called by all solvers in mesh_init in order to initialize the ip coordinates. -! Later on the current ip coordinates are directly prvided by the spectral solver and by Abaqus, -! so no need to use this subroutine anymore; Marc however only provides nodal displacements, -! so in this case the ip coordinates are always calculated on the basis of this subroutine. -! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! FOR THE MOMENT THIS SUBROUTINE ACTUALLY CALCULATES THE CELL CENTER AND NOT THE IP COORDINATES, -! AS THE IP IS NOT (ALWAYS) LOCATED IN THE CENTER OF THE IP VOLUME. -! HAS TO BE CHANGED IN A LATER VERSION. -! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -!-------------------------------------------------------------------------------------------------- -subroutine mesh_build_ipCoordinates + do e = 1,size(connectivity,3) + do i = 1,elem%nIPs - integer :: e,i,n - real(pReal), dimension(3) :: myCoords + select case (elem%cellType) + case (1) ! 2D 3node + IPvolume(i,e) = math_areaTriangle(node(1:3,connectivity(1,i,e)), & + node(1:3,connectivity(2,i,e)), & + node(1:3,connectivity(3,i,e))) - !$OMP PARALLEL DO PRIVATE(myCoords) - do e = 1,theMesh%nElems - do i = 1,theMesh%elem%nIPs - myCoords = 0.0_pReal - do n = 1,theMesh%elem%nCellnodesPerCell - myCoords = myCoords + mesh_cellnode(1:3,mesh_cell2(n,i,e)) - enddo - mesh_ipCoordinates(1:3,i,e) = myCoords / real(theMesh%elem%nCellnodesPerCell,pReal) + case (2) ! 2D 4node + IPvolume(i,e) = math_areaTriangle(node(1:3,connectivity(1,i,e)), & ! assume planar shape, division in two triangles suffices + node(1:3,connectivity(2,i,e)), & + node(1:3,connectivity(3,i,e))) & + + math_areaTriangle(node(1:3,connectivity(3,i,e)), & + node(1:3,connectivity(4,i,e)), & + node(1:3,connectivity(1,i,e))) + case (3) ! 3D 4node + IPvolume(i,e) = math_volTetrahedron(node(1:3,connectivity(1,i,e)), & + node(1:3,connectivity(2,i,e)), & + node(1:3,connectivity(3,i,e)), & + node(1:3,connectivity(4,i,e))) + case (4) ! 3D 8node + ! J. Grandy, Efficient Calculation of Volume of Hexahedral Cells + ! Lawrence Livermore National Laboratory + ! https://www.osti.gov/servlets/purl/632793 + x0 = node(1:3,connectivity(1,i,e)) + x1 = node(1:3,connectivity(2,i,e)) + x2 = node(1:3,connectivity(4,i,e)) + x3 = node(1:3,connectivity(3,i,e)) + x4 = node(1:3,connectivity(5,i,e)) + x5 = node(1:3,connectivity(6,i,e)) + x6 = node(1:3,connectivity(8,i,e)) + x7 = node(1:3,connectivity(7,i,e)) + IPvolume(i,e) = dot_product((x7-x1)+(x6-x0),math_cross((x7-x2), (x3-x0))) & + + dot_product((x6-x0), math_cross((x7-x2)+(x5-x0),(x7-x4))) & + + dot_product((x7-x1), math_cross((x5-x0), (x7-x4)+(x3-x0))) + IPvolume(i,e) = IPvolume(i,e)/12.0_pReal + end select enddo enddo - !$OMP END PARALLEL DO -end subroutine mesh_build_ipCoordinates +end function IPvolume !-------------------------------------------------------------------------------------------------- -!> @brief calculation of IP interface areas, allocate globals '_ipArea', and '_ipAreaNormal' +!> @brief calculation of IP interface areas !-------------------------------------------------------------------------------------------------- -subroutine mesh_build_ipAreas +function IPareaNormal(elem,nElem,connectivity,node) - integer :: e,t,g,c,i,f,n,m - real(pReal), dimension (3,FE_maxNcellnodesPerCellface) :: nodePos, normals - real(pReal), dimension(3) :: normal + type(tElement), intent(in) :: elem + integer, intent(in) :: nElem + integer, dimension(:,:,:), intent(in) :: connectivity + real(pReal), dimension(:,:), intent(in) :: node + + real(pReal), dimension(3,elem%nIPneighbors,elem%nIPs,nElem) :: ipAreaNormal - allocate(mesh_ipArea(theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) - allocate(mesh_ipAreaNormal(3,theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) + real(pReal), dimension (3,size(elem%cellFace,1)) :: nodePos + integer :: e,i,f,n,m + + m = size(elem%cellFace,1) - !$OMP PARALLEL DO PRIVATE(t,g,c,nodePos,normal,normals) - do e = 1,theMesh%nElems ! loop over cpElems - t = mesh_element(2,e) ! get element type - g = theMesh%elem%geomType - c = theMesh%elem%cellType - select case (c) + do e = 1,nElem + do i = 1,elem%nIPs + do f = 1,elem%nIPneighbors + nodePos = node(1:3,connectivity(elem%cellface(1:m,f),i,e)) + + select case (elem%cellType) + case (1,2) ! 2D 3 or 4 node + IPareaNormal(1,f,i,e) = nodePos(2,2) - nodePos(2,1) ! x_normal = y_connectingVector + IPareaNormal(2,f,i,e) = -(nodePos(1,2) - nodePos(1,1)) ! y_normal = -x_connectingVector + IPareaNormal(3,f,i,e) = 0.0_pReal + case (3) ! 3D 4node + IPareaNormal(1:3,f,i,e) = math_cross(nodePos(1:3,2) - nodePos(1:3,1), & + nodePos(1:3,3) - nodePos(1:3,1)) + case (4) ! 3D 8node + ! for this cell type we get the normal of the quadrilateral face as an average of + ! four normals of triangular subfaces; since the face consists only of two triangles, + ! the sum has to be divided by two; this whole prcedure tries to compensate for + ! probable non-planar cell surfaces + IPareaNormal(1:3,f,i,e) = 0.0_pReal + do n = 1, m + IPareaNormal(1:3,f,i,e) = IPareaNormal(1:3,f,i,e) & + + math_cross(nodePos(1:3,mod(n+0,m)+1) - nodePos(1:3,n), & + nodePos(1:3,mod(n+1,m)+1) - nodePos(1:3,n)) * 0.5_pReal + enddo + end select + enddo + enddo + enddo - case (1,2) ! 2D 3 or 4 node - do i = 1,theMesh%elem%nIPs - do f = 1,FE_NipNeighbors(c) ! loop over cell faces - forall(n = 1:FE_NcellnodesPerCellface(c)) & - nodePos(1:3,n) = mesh_cellnode(1:3,mesh_cell(theMesh%elem%cellface(n,f),i,e)) - normal(1) = nodePos(2,2) - nodePos(2,1) ! x_normal = y_connectingVector - normal(2) = -(nodePos(1,2) - nodePos(1,1)) ! y_normal = -x_connectingVector - normal(3) = 0.0_pReal - mesh_ipArea(f,i,e) = norm2(normal) - mesh_ipAreaNormal(1:3,f,i,e) = normal / norm2(normal) ! ensure unit length of area normal - enddo - enddo - - case (3) ! 3D 4node - do i = 1,theMesh%elem%nIPs - do f = 1,FE_NipNeighbors(c) ! loop over cell faces - forall(n = 1:FE_NcellnodesPerCellface(c)) & - nodePos(1:3,n) = mesh_cellnode(1:3,mesh_cell(theMesh%elem%cellface(n,f),i,e)) - normal = math_cross(nodePos(1:3,2) - nodePos(1:3,1), & - nodePos(1:3,3) - nodePos(1:3,1)) - mesh_ipArea(f,i,e) = norm2(normal) - mesh_ipAreaNormal(1:3,f,i,e) = normal / norm2(normal) ! ensure unit length of area normal - enddo - enddo - - case (4) ! 3D 8node - ! for this cell type we get the normal of the quadrilateral face as an average of - ! four normals of triangular subfaces; since the face consists only of two triangles, - ! the sum has to be divided by two; this whole prcedure tries to compensate for - ! probable non-planar cell surfaces - m = FE_NcellnodesPerCellface(c) - do i = 1,theMesh%elem%nIPs - do f = 1,FE_NipNeighbors(c) ! loop over cell faces - forall(n = 1:FE_NcellnodesPerCellface(c)) & - nodePos(1:3,n) = mesh_cellnode(1:3,mesh_cell(theMesh%elem%cellface(n,f),i,e)) - forall(n = 1:FE_NcellnodesPerCellface(c)) & - normals(1:3,n) = 0.5_pReal & - * math_cross(nodePos(1:3,1+mod(n ,m)) - nodePos(1:3,n), & - nodePos(1:3,1+mod(n+1,m)) - nodePos(1:3,n)) - normal = 0.5_pReal * sum(normals,2) - mesh_ipArea(f,i,e) = norm2(normal) - mesh_ipAreaNormal(1:3,f,i,e) = normal / norm2(normal) - enddo - enddo - - end select - enddo - !$OMP END PARALLEL DO - -end subroutine mesh_build_ipAreas +end function IPareaNormal !-------------------------------------------------------------------------------------------------- diff --git a/src/plastic_isotropic.f90 b/src/plastic_isotropic.f90 index 2e9a6b57e..18b6f8e86 100644 --- a/src/plastic_isotropic.f90 +++ b/src/plastic_isotropic.f90 @@ -434,7 +434,7 @@ end function plastic_isotropic_postResults !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- subroutine plastic_isotropic_results(instance,group) -#if defined(PETSc) || defined(DAMASKHDF5) +#if defined(PETSc) || defined(DAMASK_HDF5) integer, intent(in) :: instance character(len=*), intent(in) :: group diff --git a/src/plastic_nonlocal.f90 b/src/plastic_nonlocal.f90 index db5e3f097..e7e4e713d 100644 --- a/src/plastic_nonlocal.f90 +++ b/src/plastic_nonlocal.f90 @@ -1452,9 +1452,9 @@ subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & opposite_n, & !< neighbor index pointing to me when looking from my opposite neighbor t, & !< type of dislocation o,& !< offset shortcut - no,& !< neighbour offset shortcut + no,& !< neighbor offset shortcut p,& !< phase shortcut - np,& !< neighbour phase shortcut + np,& !< neighbor phase shortcut topp, & !< type of dislocation with opposite sign to t s !< index of my current slip system real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),10) :: & @@ -1653,7 +1653,7 @@ subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & !* If it's not at all compatible, no flux is arriving, because everything is dammed in front of !* my neighbor's interface. !* The entering flux from my neighbor will be distributed on my slip systems according to the - !*compatibility + !* compatibility considerEnteringFlux = .false. neighbor_v = 0.0_pReal ! needed for check of sign change in flux density below diff --git a/src/results.f90 b/src/results.f90 index d13a36fde..36b122776 100644 --- a/src/results.f90 +++ b/src/results.f90 @@ -65,7 +65,7 @@ subroutine results_init write(6,'(/,a)') ' <<<+- results init -+>>>' write(6,'(/,a)') ' Diehl et al., Integrating Materials and Manufacturing Innovation 6(1):83–91, 2017' - write(6,'(a)') ' https://doi.org/10.1007/s40192-018-0118-7' + write(6,'(a)') ' https://doi.org/10.1007/s40192-017-0084-5' resultsFile = HDF5_openFile(trim(getSolverJobName())//'.hdf5','w',.true.) call HDF5_addAttribute(resultsFile,'DADF5-version',0.3_pReal) @@ -296,21 +296,34 @@ end subroutine results_writeVectorDataset_real !-------------------------------------------------------------------------------------------------- !> @brief stores a tensor dataset in a group !-------------------------------------------------------------------------------------------------- -subroutine results_writeTensorDataset_real(group,dataset,label,description,SIunit) +subroutine results_writeTensorDataset_real(group,dataset,label,description,SIunit,transposed) character(len=*), intent(in) :: label,group,description character(len=*), intent(in), optional :: SIunit + logical, intent(in), optional :: transposed real(pReal), intent(in), dimension(:,:,:) :: dataset integer :: i + logical :: T integer(HID_T) :: groupHandle real(pReal), dimension(:,:,:), allocatable :: dataset_transposed + + if(present(transposed)) then + T = transposed + else + T = .true. + endif - allocate(dataset_transposed,mold=dataset) - do i=1,size(dataset,3) - dataset_transposed(1:3,1:3,i) = transpose(dataset(1:3,1:3,i)) - enddo + if(T) then + if(size(dataset,1) /= size(dataset,2)) call IO_error(0,ext_msg='transpose non-symmetric tensor') + allocate(dataset_transposed,mold=dataset) + do i=1,size(dataset_transposed,3) + dataset_transposed(:,:,i) = transpose(dataset(:,:,i)) + enddo + else + allocate(dataset_transposed,source=dataset) + endif groupHandle = results_openGroup(group)